Analyses for Graded encoding of spatial novelty scales in the human brain
Joern Alexander Quent
2025-10-30
- Preparation
- Results
- Methods
- Supplementary material
- Supplementary notes
- Supplementary tables
- Supplementary figures
- Learning curves
- Locomotion by novelty
- Spatial distribution of novelty scores
- Minimal influence of centrality on spatial novelty effects
- Individual differences in centrality do not predict spatial novelty
- Hippocampal gradient based on maximum statistic
- Hippocampal gradient without averaging
- Relationship between spatial memory and successful memory encoding
- Distribution of principal functional connectivity gradient scores across spatial novelty preference levels in the posterior medial cortex
- Hippocampal connectivity patterns correlating with Margulies gradient
- Posterior predictive checks for the novelty score model
- Distribuion of starting relative to object locations
- Correlation between novelty measures
Preparation
Libs and parameters
# Libraries
library(ggplot2)
# install via devtools:::install_github("JAQquent/assortedRFunctions", upgrade = "never")
library(assortedRFunctions)
library(cowplot)
library(stringr)
library(plyr)
library(mgcv)
library(foreach)
library(doParallel)
library(BayesFactor)
library(ciftiTools)
library(viridis)
library(tidybayes)
library(bayesplot)
library(brms)
library(data.table)
library(knitr)
library(readxl)
library(ggridges)
library(ggforce)
library(marginaleffects)
library(gghalves)
library(effectsize)
library(rstan)
library(latex2exp)
# For SVM
library(caTools)
library(e1071)
# Load ciftiTools and set workbench paths
possible_wb_paths <- c("/usr/bin/wb_command", "/home1/Jaquent/Toolboxes/workbench/bin_rh_linux64/")
load_ciftiTools(possible_wb_paths)Some settings that change from computer to computer that I work on.
# Use correct locations and other settings based on computer
if(Sys.info()[4] == "DESKTOP-335I26I"){
# Work laptop (Windows)
path2imaging_results2 <- "D:/Seafile/imaging_results"
} else if(Sys.info()[4] == 'DESKTOP-91CQCSQ') {
# Work desktop (Windows)
path2imaging_results2 <- "D:/imaging_results"
} else if(Sys.info()[4] == 'alex-Zenbook-UX3404VA-UX3404VA') {
# Work laptop (Linux)
path2imaging_results2 <- "/media/alex/shared/Seafile/imaging_results"
} else if(Sys.info()[4] == "greengoblin"){
# Main desktop PC (Linux)
path2imaging_results2 <- "/media/alex/work/Seafile/imaging_results"
} else if(Sys.info()[4] == "GREEN-GOBLIN-WI"){
# Main desktop PC
path2imaging_results2 <- "E:/Seafile/imaging_results"
} else {
# Personal laptop (Windows)
path2imaging_results2 <- "D:/OLM/imaging_results"
}
# Seed
set.seed(20240205)Click here for detailed session information.
sessioninfo::session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.4.0 (2024-04-24)
## os Ubuntu 22.04.5 LTS
## system x86_64, linux-gnu
## ui X11
## language en_GB:en
## collate en_GB.UTF-8
## ctype en_GB.UTF-8
## tz Asia/Shanghai
## date 2025-10-30
## pandoc 3.1.1 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## abind 1.4-5 2016-07-21 [1] CRAN (R 4.4.0)
## arrayhelpers 1.1-0 2020-02-04 [1] CRAN (R 4.4.0)
## assortedRFunctions * 1.1.3 2025-09-28 [1] Github (JAquent/assortedRFunctions@a84af63)
## backports 1.5.0 2024-05-23 [1] CRAN (R 4.4.0)
## base64enc 0.1-3 2015-07-28 [1] CRAN (R 4.4.0)
## BayesFactor * 0.9.12-4.7 2024-01-24 [1] CRAN (R 4.4.0)
## bayesplot * 1.11.1 2024-02-15 [1] CRAN (R 4.4.0)
## bayestestR 0.16.1 2025-07-01 [1] CRAN (R 4.4.0)
## bitops 1.0-7 2021-04-24 [1] CRAN (R 4.4.0)
## bridgesampling 1.1-2 2021-04-16 [1] CRAN (R 4.4.0)
## brms * 2.22.0 2024-09-23 [1] CRAN (R 4.4.0)
## Brobdingnag 1.2-9 2022-10-19 [1] CRAN (R 4.4.0)
## bslib 0.7.0 2024-03-29 [1] CRAN (R 4.4.0)
## cachem 1.1.0 2024-05-16 [1] CRAN (R 4.4.0)
## caTools * 1.18.2 2021-03-28 [1] CRAN (R 4.4.0)
## cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.4.0)
## checkmate 2.3.1 2023-12-04 [1] CRAN (R 4.4.0)
## ciftiTools * 0.17.4 2025-01-23 [1] CRAN (R 4.4.0)
## class 7.3-20 2022-01-13 [4] CRAN (R 4.1.2)
## cli 3.6.2 2023-12-11 [1] CRAN (R 4.4.0)
## coda * 0.19-4.1 2024-01-31 [1] CRAN (R 4.4.0)
## codetools 0.2-18 2020-11-04 [4] CRAN (R 4.0.3)
## colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.4.0)
## cowplot * 1.1.3 2024-01-22 [1] CRAN (R 4.4.0)
## data.table * 1.17.8 2025-07-10 [1] CRAN (R 4.4.0)
## datawizard 1.1.0 2025-05-09 [1] CRAN (R 4.4.0)
## digest 0.6.35 2024-03-11 [1] CRAN (R 4.4.0)
## distributional 0.4.0 2024-02-07 [1] CRAN (R 4.4.0)
## doParallel * 1.0.17 2022-02-07 [1] CRAN (R 4.4.0)
## dplyr 1.1.4 2023-11-17 [1] CRAN (R 4.4.0)
## e1071 * 1.7-14 2023-12-06 [1] CRAN (R 4.4.0)
## effectsize * 0.8.8 2024-05-12 [1] CRAN (R 4.4.0)
## emmeans 1.10.4 2024-08-21 [1] CRAN (R 4.4.0)
## estimability 1.5.1 2024-05-12 [1] CRAN (R 4.4.0)
## evaluate 0.23 2023-11-01 [1] CRAN (R 4.4.0)
## fansi 1.0.6 2023-12-08 [1] CRAN (R 4.4.0)
## farver 2.1.2 2024-05-13 [1] CRAN (R 4.4.0)
## fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.4.0)
## foreach * 1.5.2 2022-02-02 [1] CRAN (R 4.4.0)
## generics 0.1.3 2022-07-05 [1] CRAN (R 4.4.0)
## ggdist 3.3.2 2024-03-05 [1] CRAN (R 4.4.0)
## ggforce * 0.4.2 2024-02-19 [1] CRAN (R 4.4.0)
## gghalves * 0.1.4 2022-11-20 [1] CRAN (R 4.4.0)
## ggplot2 * 3.5.2 2025-04-09 [1] CRAN (R 4.4.0)
## ggridges * 0.5.6 2024-01-23 [1] CRAN (R 4.4.0)
## gifti 0.8.0 2020-11-11 [1] CRAN (R 4.4.0)
## glue 1.7.0 2024-01-09 [1] CRAN (R 4.4.0)
## gridExtra 2.3 2017-09-09 [1] CRAN (R 4.4.0)
## gtable 0.3.5 2024-04-22 [1] CRAN (R 4.4.0)
## htmltools 0.5.8.1 2024-04-04 [1] CRAN (R 4.4.0)
## inline 0.3.19 2021-05-31 [1] CRAN (R 4.4.0)
## insight 1.3.1 2025-06-30 [1] CRAN (R 4.4.0)
## iterators * 1.0.14 2022-02-05 [1] CRAN (R 4.4.0)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.4.0)
## jsonlite 1.8.8 2023-12-04 [1] CRAN (R 4.4.0)
## knitr * 1.50 2025-03-16 [1] CRAN (R 4.4.0)
## latex2exp * 0.9.6 2022-11-28 [1] CRAN (R 4.4.0)
## lattice 0.20-45 2021-09-22 [4] CRAN (R 4.1.1)
## lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.4.0)
## loo 2.8.0 2024-07-03 [1] CRAN (R 4.4.0)
## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.4.0)
## marginaleffects * 0.28.0 2025-06-25 [1] CRAN (R 4.4.0)
## MASS 7.3-64 2025-01-04 [1] CRAN (R 4.4.0)
## Matrix * 1.7-0 2024-04-26 [1] CRAN (R 4.4.0)
## MatrixModels 0.5-3 2023-11-06 [1] CRAN (R 4.4.0)
## matrixStats 1.3.0 2024-04-11 [1] CRAN (R 4.4.0)
## mgcv * 1.8-39 2022-02-24 [4] CRAN (R 4.1.2)
## multcomp 1.4-26 2024-07-18 [1] CRAN (R 4.4.0)
## munsell 0.5.1 2024-04-01 [1] CRAN (R 4.4.0)
## mvtnorm 1.2-5 2024-05-21 [1] CRAN (R 4.4.0)
## nlme * 3.1-155 2022-01-13 [4] CRAN (R 4.1.2)
## oro.nifti 0.11.4 2022-08-10 [1] CRAN (R 4.4.0)
## parameters 0.27.0 2025-07-09 [1] CRAN (R 4.4.0)
## pbapply 1.7-2 2023-06-27 [1] CRAN (R 4.4.0)
## pillar 1.9.0 2023-03-22 [1] CRAN (R 4.4.0)
## pkgbuild 1.4.4 2024-03-17 [1] CRAN (R 4.4.0)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.4.0)
## plyr * 1.8.9 2023-10-02 [1] CRAN (R 4.4.0)
## polyclip 1.10-6 2023-09-27 [1] CRAN (R 4.4.0)
## posterior 1.6.1 2025-02-27 [1] CRAN (R 4.4.0)
## proxy 0.4-27 2022-06-09 [1] CRAN (R 4.4.0)
## purrr 1.0.2 2023-08-10 [1] CRAN (R 4.4.0)
## QuickJSR 1.2.2 2024-06-07 [1] CRAN (R 4.4.0)
## R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.4.0)
## R.oo 1.26.0 2024-01-24 [1] CRAN (R 4.4.0)
## R.utils 2.12.3 2023-11-18 [1] CRAN (R 4.4.0)
## R6 2.5.1 2021-08-19 [1] CRAN (R 4.4.0)
## RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.4.0)
## Rcpp * 1.0.12 2024-01-09 [1] CRAN (R 4.4.0)
## RcppParallel 5.1.7 2023-02-27 [1] CRAN (R 4.4.0)
## readxl * 1.4.3 2023-07-06 [1] CRAN (R 4.4.0)
## rlang 1.1.4 2024-06-04 [1] CRAN (R 4.4.0)
## rmarkdown 2.27 2024-05-17 [1] CRAN (R 4.4.0)
## RNifti 1.6.1 2024-03-07 [1] CRAN (R 4.4.0)
## rstan * 2.32.6 2024-03-05 [1] CRAN (R 4.4.0)
## rstantools 2.4.0 2024-01-31 [1] CRAN (R 4.4.0)
## rstudioapi 0.16.0 2024-03-24 [1] CRAN (R 4.4.0)
## sandwich 3.1-1 2024-09-15 [1] CRAN (R 4.4.0)
## sass 0.4.9 2024-03-15 [1] CRAN (R 4.4.0)
## scales 1.3.0 2023-11-28 [1] CRAN (R 4.4.0)
## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.4.0)
## StanHeaders * 2.32.9 2024-05-29 [1] CRAN (R 4.4.0)
## stringi 1.8.4 2024-05-06 [1] CRAN (R 4.4.0)
## stringr * 1.5.1 2023-11-14 [1] CRAN (R 4.4.0)
## survival 3.2-13 2021-08-24 [4] CRAN (R 4.1.1)
## svUnit 1.0.6 2021-04-19 [1] CRAN (R 4.4.0)
## tensorA 0.36.2.1 2023-12-13 [1] CRAN (R 4.4.0)
## TH.data 1.1-2 2023-04-17 [1] CRAN (R 4.4.0)
## tibble 3.2.1 2023-03-20 [1] CRAN (R 4.4.0)
## tidybayes * 3.0.6 2023-08-12 [1] CRAN (R 4.4.0)
## tidyr 1.3.1 2024-01-24 [1] CRAN (R 4.4.0)
## tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.4.0)
## tweenr 2.0.3 2024-02-26 [1] CRAN (R 4.4.0)
## utf8 1.2.4 2023-10-22 [1] CRAN (R 4.4.0)
## vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.4.0)
## viridis * 0.6.5 2024-01-29 [1] CRAN (R 4.4.0)
## viridisLite * 0.4.2 2023-05-02 [1] CRAN (R 4.4.0)
## withr 3.0.0 2024-01-16 [1] CRAN (R 4.4.0)
## xfun 0.52 2025-04-02 [1] CRAN (R 4.4.0)
## xml2 1.3.6 2023-12-04 [1] CRAN (R 4.4.0)
## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.4.0)
## yaml 2.3.8 2023-12-11 [1] CRAN (R 4.4.0)
## zoo 1.8-12 2023-04-13 [1] CRAN (R 4.4.0)
##
## [1] /home/alex/R/x86_64-pc-linux-gnu-library/4.4
## [2] /usr/local/lib/R/site-library
## [3] /usr/lib/R/site-library
## [4] /usr/lib/R/library
##
## ──────────────────────────────────────────────────────────────────────────────Click here for chunk for figure and other reporting parameters.
# Significance cut off
cutOff <- 1.301 # equals -log(0.05, 10)
# Parameters how to report means
report_type <- 1
digits1 <- 2
rounding_type <- "signif"
theme_journal <- function(base_size = 7, linewidth = 0.35, theme_name = "classic") {
# Select a base theme from https://ggplot2.tidyverse.org/reference/ggtheme.html
if(theme_name == "grey" | theme_name == "gray"){
base_theme <- theme_grey(base_size = base_size)
} else if(theme_name == "bw"){
base_theme <- theme_bw(base_size = base_size)
} else if(theme_name == "linedraw"){
base_theme <- theme_linedraw(base_size = base_size)
} else if(theme_name == "light"){
base_theme <- theme_light(base_size = base_size)
} else if(theme_name == "dark"){
base_theme <- theme_dark(base_size = base_size)
} else if(theme_name == "minimal"){
base_theme <- theme_minimal(base_size = base_size)
} else if(theme_name == "classic"){
base_theme <- theme_classic(base_size = base_size)
} else {
stop("Unknown theme. Check https://ggplot2.tidyverse.org/reference/ggtheme.html")
}
# Base them plus customisation
base_theme +
theme(text = element_text(color = "black", size = base_size),
axis.text = element_text(size = base_size * 0.9), # 90% of base
axis.line = element_line(linewidth = linewidth),
axis.ticks = element_line(linewidth = linewidth),
legend.text = element_text(size = base_size * 0.9), # 90% of base
legend.key.size = unit(0.5, "lines"),
plot.title = element_text(hjust = 0.5, size = base_size))
}
# Updating geom's defaults
update_geom_defaults("point", list(size = 0.5, stroke = 0.25))
update_geom_defaults("line", list(linewidth = 0.10))
update_geom_defaults("density", list(linewidth = 0.25))
update_geom_defaults("boxplot", list(linewidth = 0.25, size = 0.2))
update_geom_defaults("smooth", list(linewidth = 0.5))
update_geom_defaults("hline", list(linewidth = 0.25))
update_geom_defaults("vline", list(linewidth = 0.25))
update_geom_defaults("curve", list(linewidth = 0.25))
update_geom_defaults("segment", list(linewidth = 0.25))
# Parameters for saving figures
## Information: https://www.nature.com/ncomms/submit/how-to-submit
### PDF page: 210 x 276 mm
### Column widths in mm
single_column <- 88
double_column <- 180
dpi <- 1000
figurePath <- "figures/SpaNov/"
# Colours used for visualisation
boxplot_border_col <- "black"
boxplot_point_col <- "darkgrey"
boxplot_line_col <- "darkgrey"
base_col <- c("#7091C0", "#4A66AC", "#364875")
spectral_colours <- c("#5956a5", "#a60a44")
cool_warm_colours <- c("#3C4DC1", "#B70B28")
novFam_gradient <- viridis(n = 6, option = "H", direction = -1)
# Information for Yeo 7
## Names etc. for the networks in Yeo 7
Yeo7_fullNames <- c("Frontoparietal", "Default", "Dorsal Attention", "Limbic",
"Ventral Attention", "Somatomotor", "Visual")
Yeo7_abbr <- c("Cont", "Default", "DorsAttn", "Limbic", "SalVentAttn", "SomMot", "Vis")
## Colours used for Yeo 7
Yeo7_colours <- c("#E79523", "#CD3E4E", "#00760F", "#DCF8A4", "#C43BFA", "#4682B4", "#781286")Cifti files
Click here for support CIFTI files etc.
# Loading the support CIFTI files
## Place where to find some of the CIFTI files and parcellations
CIFTI_locations <- "data/ignore_fMRI_version1/sourceFiles/"
## CIFTI files
parcellationFile <- "Q1-Q6_RelatedValidation210.CorticalAreas_dil_Final_Final_Areas_Group_Colors_with_Atlas_ROIs2.32k_fs_LR.dlabel.nii"
CAB_NP <- "CortexSubcortex_ColeAnticevic_NetPartition_wSubcorGSR_netassignments_LR.dlabel.nii"
surfLeft <- "S1200.L.inflated_MSMAll.32k_fs_LR.surf.gii"
surfRight <- "S1200.R.inflated_MSMAll.32k_fs_LR.surf.gii"
## Combine CIFTI_locations with file name
parcellationFile <- paste0(CIFTI_locations, parcellationFile)
CAB_NP <- paste0(CIFTI_locations, CAB_NP)
surfLeft <- paste0(CIFTI_locations, surfLeft)
surfRight <- paste0(CIFTI_locations, surfRight)
## Loading CIFTIw via ciftiTools as xiis
### Get MMP parcellation
MMP_xii <- ciftiTools::read_cifti(parcellationFile,
brainstructures = "all",
surfL_fname = surfLeft,
surfR_fname = surfRight)
### Load Yeo 7 parcellation
Yeo7_xii <- ciftiTools::read_cifti("other_stuff/Yeo7.dlabel.nii",
surfL_fname = surfLeft,
surfR_fname = surfRight)
## Load other stuff
### Load the parcel names for MMP
parcel_names <- read.csv("data/ignore_fMRI_version1/extracted_values/Parcellations/MP1.0_210V_parcel_names.csv", header = FALSE)
### Load the extracted MNI coordinates
MNI_coord <- read.csv("other_stuff/cifti_subcortical_MNI152_coordinates.csv")
### Load the hippocampal projection values
load("other_stuff/projected_HC.RData")Custom functions
Click here for code for custom functions.
create_str_from_avg_comparisons <- function (x, measure_name = "X" , digits = 2, rounding_type = "round"){
# Convert to numeric values
x_num <- suppressWarnings(as.numeric(x))
# Round and create string
if (rounding_type == "signif") {
return_string <- paste0(signif(x_num[3], digits), " ", measure_name , " (95 % CI [",
signif(x_num[4], digits), ", ", signif(x_num[5], digits),
"])")
}
else if (rounding_type == "round") {
return_string <- paste0(round(x_num[3], digits), " ", measure_name , " (95 % CI [",
round(x_num[4], digits), ", ", round(x_num[5], digits),
"])")
}
else {
stop("Wrong rounding type. Choose signif or round.")
}
return(return_string)
}Loading and preparing the data
Click here for loading data.
# Load the data
## Specify paths where the data is saved
path2data <- "data/ignore_fMRI_version1/"
EV_folder <- "ignore_eventTable2/"
## Load the look-up table that contains information of R-numbers which are retracted
lookupTable <- read.csv(paste0(path2data, "lookUpTable.csv"))
## Load .RData images of the combined data (all subject in one DF)
load(paste0(path2data, "combined_data/demographics.RData"))
load(paste0(path2data, "combined_data/DW_all_data.RData"))
load(paste0(path2data, "combined_data/OLM_7T_all_data.RData"))
load(paste0(path2data, "combined_data/OLM_3T_all_data.RData"))
load(paste0(path2data, "combined_data/question_data.RData"))
# Select the subjects included in this analysis
## Load the subjects that are included in this analysis
subjectFile <- readLines(paste0(path2data, "SpaNov_subject2analyse.txt"))
subjIDs_R <- str_split(subjectFile, pattern = ",")[[1]]
subjIDs <- lookupTable$anonKey[lookupTable$Rnum %in% subjIDs_R]
# Important note: subjIDs_R do not have the same order as subjIDs!!!!!!!!!!!!!
## Subset to data that is being included in the analysis
OLM_7T_position_data <- OLM_7T_position_data[OLM_7T_position_data$subject %in% subjIDs, ]
demographics <- demographics[demographics$subject %in% subjIDs, ]
DW_position_data <- DW_position_data[DW_position_data$subject %in% subjIDs, ]
OLM_7T_logEntries <- OLM_7T_logEntries[OLM_7T_logEntries$ppid %in% subjIDs, ]
OLM_7T_trial_results <- OLM_7T_trial_results[OLM_7T_trial_results$subject %in% subjIDs, ]
question_data <- question_data[question_data$subject %in% subjIDs, ]
# Create positionData
positionData <- OLM_7T_position_data
# Add ppid to this data frame
subjects_anon <- unique(positionData$subject)
positionData$ppid <- NA
for(i in 1:length(subjects_anon)){
positionData$ppid[positionData$subject == subjects_anon[i]] <- lookupTable$Rnum[lookupTable$anonKey == subjects_anon[i]]
}
# Add ppid to OLM_7T_trial_results
OLM_7T_trial_results$ppid <- NA
for(i in 1:nrow(OLM_7T_trial_results)){
OLM_7T_trial_results$ppid[i] <- lookupTable$Rnum[lookupTable$anonKey == OLM_7T_trial_results$subject[i]]
}
# Subset to retrieval only
OLM_7T_retrieval <- OLM_7T_trial_results[OLM_7T_trial_results$trialType == "retrieval", ]
# Get the object locations to verify object placement in screenshots
obj_locations <- ddply(OLM_7T_trial_results, c("targets", "objectName", "targetNames"),
summarise, object_x_sd = sd(object_x), object_x = mean(object_x),
object_z_sd = sd(object_z), object_z = mean(object_z))
# Subset to encoding only
OLM_7T_encoding <- OLM_7T_trial_results[OLM_7T_trial_results$trialType == "encoding", ]
# Get cue times
cues_agg <- data.frame(ppid = OLM_7T_encoding$ppid,
trial = OLM_7T_encoding$trial_num,
start = OLM_7T_encoding$start_time,
end = OLM_7T_encoding$start_time + OLM_7T_encoding$cue,
duration = OLM_7T_encoding$cue,
block_num = OLM_7T_encoding$block_num)
# Get delay times
delays_agg <- data.frame(ppid = OLM_7T_encoding$ppid,
trial = OLM_7T_encoding$trial_num,
start = OLM_7T_encoding$start_time + OLM_7T_encoding$cue,
end = OLM_7T_encoding$start_time + OLM_7T_encoding$cue + OLM_7T_encoding$delay,
duration = OLM_7T_encoding$delay,
block_num = OLM_7T_encoding$block_num)# Parameter
downsample_value <- 0.2 #200 msec/0.2 sec
# Add run information
run_info <- ddply(OLM_7T_trial_results, c("block_num", "trial_num"), summarise, N = length(ppid))
positionData$run <- NA
# Loop through all trials
for(i in 1:nrow(run_info)){
positionData$run[positionData$trial == run_info$trial_num[i]] <- run_info$block_num[i]
}
# Prepare the cluster for a parallel loop
# Create the cluster following https://www.blasbenito.com/post/02_parallelizing_loops_with_r/
my.cluster <- parallel::makeCluster(detectCores() - 2, type = "PSOCK")
# Register it to be used by %dopar%
doParallel::registerDoParallel(cl = my.cluster)
# Algorithm to down sample to a sample every x msec
# For that loop through tempData_currentSubj and check if 200 msec have passed since current time
# Each row where a new run begins is always included. After that only a sample that is 200 msec
# passed the current time. Especially the first sample of a run is important to accurately
# determine when the first image is recorded and that all the onsets correspond to that.
# Extract unique participants but in order in which they appear. This was verified
# E.g. via positionData$ppid[!duplicated(positionData$ppid)]
subjects <- unique(positionData$ppid)
# Run parallel loop
include <- foreach(i = 1:length(subjIDs_R), .combine = 'c') %dopar% {
# Subset to current subject
tempData_currentSubj <- positionData[positionData$ppid == subjIDs_R[i], ]
# Create include variable
tempInclude <- rep(FALSE, nrow(tempData_currentSubj))
# Loop through the each row
for(i in 1:nrow(tempData_currentSubj)){
# First time
if(i == 1){
currentTime <- tempData_currentSubj$time[i]
tempInclude[i] <- TRUE
} else {
# Check if the last row was from a different run
if(tempData_currentSubj$run[i - 1] != tempData_currentSubj$run[i]){
# Set new time in this case
currentTime <- tempData_currentSubj$time[i]
tempInclude[i] <- TRUE
} else {
# Check ih ith time is larger than currentTime + downsample_value
if(tempData_currentSubj$time[i] > currentTime + downsample_value){
# Set new time in this case
currentTime <- tempData_currentSubj$time[i]
tempInclude[i] <- TRUE
}
}
}
}
# Return
tempInclude
}
# Stop cluster again
parallel::stopCluster(cl = my.cluster)
# Apply downsampling to positionData
positionData2 <- positionData[include,]
# Check if there is no problem
time_diff <- ddply(positionData2, c("ppid", "run", "trial"), summarise,
mean_time_diff = mean(diff(time)),
sd_time_diff = sd(diff(time)))
# Now I manually checked the run start time of the first two participants
# Get the number of seeds of this analysis
numSeeds <- 10
limValues <- c(-90,90)
# Bin the environment
x <- positionData2$pos_x
z <- positionData2$pos_z
env_sectors <- voronoi_tessellation_grid_binning_2d(x, z,limValues, numSeeds, "hexagon",
useParallelisation = TRUE)
# Add result back to the df (sector2 is only for plotting)
positionData2$sector <- env_sectors
positionData2$sector2 <- factor(env_sectors, levels = sample(unique(env_sectors)))
# Save in intermediate data
save(positionData2, file = "intermediate_data/positionData2.RData")# Load data
load("intermediate_data/positionData2.RData")
# Exclude control trials
no_control <- positionData2[positionData2$trialType != "control", ]
# Include only Run 1, 2 and 3 because the last retrieval run
# doesn't count as far as this analysis is concerned
no_control <- no_control[no_control$run %in% 1:3, ]Results
Modelling of spatial novelty during naturalistic navigation
Reporting demographic information
n <- nrow(demographics)
str1 <- mean_SD_str2(demographics$age, type = 1, digits = digits1, rounding_type = rounding_type, measure = "years")
females <- table(demographics$gender)[1]
males <- table(demographics$gender)[2]- n = 56
- Age = M = 24 years (SD = 2.8 years)
- Gender ratio = 35/21
- Age range 19.7, 37.3
Figure 1
3D histogramms
get_seeds_from_voronoi_tessellation_hexagon <- function(limValues, numSeeds){
# Create bins like it is done in voronoi_tessellation_grid_binning_2d
# This is code directly from that function to get the coordinates of each sector
xLim <- sort(limValues)
yLim <- sort(limValues)
byValue <- ((xLim[2] - xLim[1])/(numSeeds - 1))
xRange <- seq(from = xLim[1], to = xLim[2], by = byValue)
xStepSize <- xRange[2] - xRange[1]
x_shiftValue <- xStepSize/4
yRange <- seq(from = yLim[1], to = yLim[2], by = xStepSize * (sqrt(3)/2))
yRange <- yRange[1:length(xRange)]
yRange <- yRange + (yLim[2] - max(yRange))/2
seeds <- data.frame(x = xRange[1], y = yRange)
shiftIndex <- rep(c(1, -1), length.out = nrow(seeds))
seeds$x <- seeds$x + shiftIndex * x_shiftValue
for(i in 2:length(xRange)){
seeds <- rbind(seeds, data.frame(x = xRange[i] + shiftIndex * x_shiftValue, y = yRange))
}
seeds$sector <- row.names(seeds)
return(seeds)
}This creates the data necessary to create the 3D histogram using blender.
# Load data from
load("E:/research_projects/OLM_project/analysis/ignore_eventTable3/images/SpaNov_event_file_contPM.RData")
# Convert list to data frame
data <- rbindlist(subj_list, idcol = "subject")
# Subject to character
data$subject <- as.character(data$subject)
# Subset to run 1, 2, 3 to include the first retrieval run
data_sub2 <- data[data$run != 4, ]
# Calculate average/sum per subject
sector_average <- ddply(data_sub2, c("subject","sector"), summarise,
visits = sum(visits, na.rm = TRUE),
lastVisit = mean(lastVisit, na.rm = TRUE))
# Calculate average per sector
sector_average <- ddply(sector_average, c("sector"), summarise,
visits = mean(visits, na.rm = TRUE),
lastVisit = mean(lastVisit, na.rm = TRUE))
# Bin the environment
limValues <- c(-90,90)
numSeeds <- 10 # Number of values per axis
# Create bins like it is done in voronoi_tessellation_grid_binning_2d
seeds <- get_seeds_from_voronoi_tessellation_hexagon(limValues, numSeeds)
# Add x & y coordinates to sector_visits. For this just loop the df
sector_average$x <- NA
sector_average$y <- NA
for(i in 1:nrow(sector_average)){
# Get current sector
currentSector <- sector_average$sector[i]
# Get corresponding row index in seeds
rowID <- which(seeds$sector == currentSector)
# Assign the correct coordinates based on rowID
sector_average$x[i] <- seeds$x[rowID]
sector_average$y[i] <- seeds$y[rowID]
}
# Prepare data for 3D plot for the average number of visits per persons to each sector
temp_data <- sector_average[,-3]
# Write .csv file
write.csv(x = temp_data, file = "other_stuff/source_data_files/blenderData_avg_visits.csv",
quote = FALSE, row.names = FALSE)
# Prepare data for 3D plot for the average number of visits per persons to each sector
## Remove NA values & remove wrong column
temp_data <- na.omit(sector_average)
temp_data <- temp_data[,-2]
# Write .csv file
write.csv(x = temp_data, file = "other_stuff/source_data_files/blenderData_avg_time_since_last_visit.csv",
quote = FALSE, row.names = FALSE)
# Colours
barColours <- viridisLite::viridis(n = 5, option = "D")Range of values is 75 to 111321.
Distribution of total path lengths
# Subset to only data from the encoding part
pathData <- positionData[positionData$trialType == 'encoding', ]
# Calculate the how much the perfect participant would need to travel
## Subset to only encoding trials
boolIndex <- OLM_7T_trial_results$trialType == "encoding"
OLM_7T_encoding <- OLM_7T_trial_results[boolIndex, ]
# Function to calculate the path lengths
calculate_path_length <- function(pos_x, pos_z){
# Calculate difference between points and square them
diff_X <- diff(pos_x)^2
diff_z <- diff(pos_z)^2
# Calculate the sum and take the square root
dists_travelled <- sqrt(diff_X + diff_z)
# Sum up all distances travelled to one path_length
return(sum(dists_travelled))
}
# Calculate the path traveled for each trial and each subject
# Use pathData because it only includes encoding data
path_lenghts <- ddply(pathData, c("ppid", "trial"),
summarise,
distance = calculate_path_length(pos_x, pos_z))
# Calculate the total path lengths for the whole experiment
total_path_lengths <- ddply(path_lenghts, c("ppid"),
summarise, distance = sum(distance))
# Write .csv file
write.csv(x = total_path_lengths, file = "other_stuff/source_data_files/total_path_lengths.csv",
quote = FALSE, row.names = FALSE)
# Create a box plot
total_path_lengths_plot <- ggplot(total_path_lengths, aes(x = distance)) +
geom_density(fill = base_col[1]) +
geom_jitter(aes(y = -0.001), height = 0.0005, width = 0, colour = boxplot_point_col) +
labs(y = "Density", x = "Total path length (vm)", title = "") +
scale_x_continuous(breaks = c(3000, 3300, 3600)) +
coord_cartesian(xlim = c(3000, 3600), expand = FALSE, ylim = c(-0.001*2, 0.0065)) +
theme_journal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.y = element_blank(),
plot.margin = unit(c(1, 2.7, 1, 2.7), "mm"))Distributions of event durations
# Load the event files
event_files <- list.files(path = "event_tables/OLMe_7T_SpaNov_gradient_6lvl/", recursive = TRUE, full.names = TRUE)
# Remove "per_tra.txt" and "dur_tra.txt" from list
event_files <- event_files[!str_detect(event_files, "tra.txt")]
# Remove cue and delay
event_files <- event_files[!str_detect(event_files, "cue.txt") & !str_detect(event_files, "delay.txt")]
# Create list
event_list <- list()
# Load the event files
for(i in 1:length(event_files)){
# Load current event file
temp_event_file <- read.table(event_files[i], header = FALSE, sep = "\t")
temp_event_file$V3 <- NULL # Remove the third column
names(temp_event_file) <- c("onset", "duration")
# Split the file path to get information
event_file_path_split <- str_split_1(event_files[i], pattern = "/")
temp_event_file$subject <- event_file_path_split[4]
temp_event_file$run <- event_file_path_split[5]
# To extract the level of the event extract "lvl" plus an integer from string
temp_event_file$lvl <- str_extract(event_file_path_split[8], "lvl[0-9]+")
# Add to list
event_list[[i]] <- temp_event_file
}
# Combine to data frame
event_df <- as.data.frame(rbindlist(event_list))
# Calculate offset
event_df$offset <- event_df$onset + event_df$duration
# Convert lvl label to number
event_df$lvl_num <- as.numeric(str_extract(event_df$lvl, "[0-9]+"))
# Change run label
event_df$run <- ifelse(event_df$run == "run-01", "Run 1", "Run 2")
# Write .csv file
write.csv(x = event_df, file = "other_stuff/source_data_files/event_duration.csv",
quote = FALSE, row.names = FALSE)
# Create a box plot
event_duration_plot <- ggplot(event_df, aes(x = duration)) +
geom_density(fill = base_col[1]) +
geom_jitter(aes(y = -0.1), height = 0.05, width = 0, colour = boxplot_point_col, alpha = 0.1) +
labs(y = "Density", x = "Event duration (s)", title = "") +
scale_x_continuous(breaks = c(0, 16, 32)) +
coord_cartesian(xlim = c(0, 32), expand = FALSE, ylim = c(-0.1*2, 0.8)) +
theme_journal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.y = element_blank(),
plot.margin = unit(c(1, 2.5, 1, 2.5), "mm"))Distributions of Rayleigh vector length
Calculate Rayleigh Test of Uniformity for each sector
#install.packages("circular") # If not already installed
library(circular)
# Function to calculate the test and apply it to the data frame
calc_rayleigh_test_FUN <- function(angles){
# Create circular object
circ_data <- circular(angles, units = "degrees", template = "none")
# Calculate Rayleigh test which includes R
rayleigh_test <- rayleigh.test(circ_data)
# Return
return(c(rayleigh_test$statistic, rayleigh_test$p.value))
}
# Apply test to all sectors
sector_rayleigh <- ddply(no_control, c("sector"), summarise,
resultant_length = calc_rayleigh_test_FUN(rot_y)[1],
p.value = calc_rayleigh_test_FUN(rot_y)[2])
# Write .csv file
write.csv(x = sector_rayleigh, file = "other_stuff/source_data_files/sector_rayleigh.csv",
quote = FALSE, row.names = FALSE)
sector_rayleigh_plot <- ggplot(sector_rayleigh, aes(x = resultant_length)) +
geom_histogram(colour = base_col[1], fill = base_col[1]) +
geom_vline(xintercept = mean(sector_rayleigh$resultant_length), colour = "black", linetype = 2, size = 0.5) +
scale_x_continuous(breaks = c(0, 0.5, 1)) +
scale_y_continuous(breaks = c(0, 5, 10)) +
coord_cartesian(xlim = c(-0.02, 1), ylim = c(0, 10),
expand = FALSE) +
labs(x = "Length of Rayleigh vector", y = "Count") +
theme_journal() +
theme(plot.margin = unit(c(1, 2.5, 1, 2.5), "mm"))
# Report mean
str1 <- mean_SD_str2(sector_rayleigh$resultant_length, report_type, digits1, rounding_type)Time spend in each locomotion state
Calculate overall time spent on translation, rotation and being stationary:
# Remove cue & delay periods
locomotion_data <- pathData
locomotion_data$include <- TRUE
subjs <- unique(locomotion_data$ppid)
trials <- unique(locomotion_data$trial)
for(i in 1:length(subjs)){
for(j in 1:length(trials)){
# Get cue start & delay end time
cue_start <- cues_agg$start[cues_agg$ppid == subjs[i] & cues_agg$trial == trials[j]]
delay_end <- delays_agg$end[delays_agg$ppid == subjs[i] & delays_agg$trial == trials[j]]
# Set include to false for those times
bool_index <- locomotion_data$ppid == subjs[i] &
locomotion_data$trial == trials[j] &
locomotion_data$time >= cue_start &
locomotion_data$time <= delay_end
locomotion_data$include[bool_index] <- FALSE
}
}
# Remove those times
locomotion_data <- locomotion_data[locomotion_data$include, ]
# The amount we found the values to avoid false positive
rotation_round <- 2 # round rotation values to this decimal point
# Function to determine which state a time point belongs to
what_state <- function(rot_y, moving2){
# Get angles
angle1 <- rot_y[2:length(rot_y)]
angle2 <- rot_y[1:(length(rot_y)-1)]
# Calculate the amount was rotated between the time points and then rotate
rotated <- c(NA, round(angularDifference(angle1, angle2), rotation_round))
# If rotation is zero called it stationary, otherwise rotation
tra_rot_sta <- ifelse(abs(rotated) == 0 | is.na(rotated),
'stationary', 'rotation')
# Set time point to translation based the information saved by unity
tra_rot_sta[moving2] <- 'translation'
# Return
return(tra_rot_sta)
}
# Calculate the state for each time points for each subject and each trial
locomotion_data <- ddply(locomotion_data, c("ppid", "trial"),
mutate, locomotion = what_state(rot_y, moving))
# Calculate the total percentage for each subject
locomotion_per <- ddply(locomotion_data, c("ppid"), summarise,
translation = mean(locomotion == "translation"),
rotation = mean(locomotion == "rotation"),
stationary = mean(locomotion == "stationary"))
# Convert from wide to long format
locomotion_per_long <- reshape2::melt(locomotion_per, id.vars = c("ppid"))
# Convert proportions to percentage
locomotion_per_long$value <- locomotion_per_long$value * 100
# Write .csv
write.csv(x = locomotion_per_long, file = "other_stuff/source_data_files/locomotion_per_long.csv",
quote = FALSE, row.names = FALSE)
# Plot showing percentage of locomotion
locomotion_states_plot <- ggplot(locomotion_per_long, aes(x = variable, y = value)) +
geom_line(aes(group = ppid), colour = boxplot_line_col) +
geom_point(colour = boxplot_point_col) +
geom_boxplot(colour = boxplot_border_col, outlier.shape = NA, width = 0.4,
fill = base_col[2]) +
stat_summary(geom = "point", fun = "mean", col = 'white', shape = 24,
fill = base_col[2], position = position_dodge(width = 0.75), size = 0.8) +
theme_journal() +
#scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
scale_y_continuous(breaks = c(0, 40, 80)) +
coord_cartesian(xlim = c(0.5, 3.5),
ylim = c(0, 80), expand = FALSE) +
theme(legend.position = "none", plot.margin = unit(c(1, 1, 1, 1), "mm")) +
labs(title = "", x = "Locomotion", y = "Percentage")Novelty score over time model
# Load the models
load("fitted_brms_models/SpaNov_m_noveltyScore.Rdata")
# Background: https://www.andrewheiss.com/blog/2022/11/29/conditional-marginal-marginaleffects/
# Set range
x1_range <- range(m_noveltyScore_run1$data$s_onset_rel)
x1_points <- seq(from = x1_range[1], to = x1_range[2], length.out = 3)
x2_range <- range(m_noveltyScore_run2$data$s_onset_rel)
x2_points <- seq(from = x2_range[1], to = x2_range[2], length.out = 3)
# Get subject-level predictions
pred1_df <- predictions(m_noveltyScore_run1, newdata = datagrid(s_onset_rel = x1_points, subject = unique),
by = c("s_onset_rel", "subject"))
pred2_df <- predictions(m_noveltyScore_run2, newdata = datagrid(s_onset_rel = x2_points, subject = unique),
by = c("s_onset_rel", "subject"))
pred1_df <- as.data.frame(pred1_df)
pred2_df <- as.data.frame(pred2_df)
# Calculate the minimum for each subject, which is used for colouring
pred1_df <- ddply(pred1_df, c("subject"), mutate, min = min(estimate))
pred2_df <- ddply(pred2_df, c("subject"), mutate, min = min(estimate))
# Create plots
## Run 1
p1 <- ggplot(data = pred1_df, aes(x = s_onset_rel, y = estimate, group = subject, colour = min)) +
geom_line() +
scale_color_viridis_c() + theme_journal() +
theme(legend.position = "none",
plot.margin = unit(c(1, 2, 1, 2), "mm")) +
labs(title = "Run 1", x = "Time", y = "Novelty score") +
scale_x_continuous(breaks = x1_range, labels = c("Start", "End")) +
scale_y_continuous(breaks = c(-2.0, -0.7, 0.6)) +
coord_cartesian(xlim = x1_range, ylim = c(-2.03, 0.6),
expand = FALSE)
## Run 2
p2 <- ggplot(data = pred2_df, aes(x = s_onset_rel, y = estimate, group = subject, colour = min)) +
geom_line() +
scale_color_viridis_c() + theme_journal() +
theme(legend.position = "none",
plot.margin = unit(c(1, 2, 1, 2), "mm")) +
labs(title = "Run 2", x = "Time", y = "Novelty score") +
scale_x_continuous(breaks = x2_range, labels = c("Start", "End")) +
scale_y_continuous(breaks = c(-1.5, -0.5, 0.5)) +
coord_cartesian(xlim = x2_range, ylim = c(-1.53, 0.5),
expand = FALSE)
# Combine and save
novelty_score_model_plot <- plot_grid(p1, p2, ncol = 2)Combine to one figure
# Combine
part1 <- plot_grid(total_path_lengths_plot, event_duration_plot,
sector_rayleigh_plot, locomotion_states_plot, align = "h", nrow = 1)
part2 <- plot_grid(NULL, novelty_score_model_plot)
Figure1_combined <- plot_grid(part1, NULL, part2, align = "v", nrow = 3, rel_heights = c(1, 0.2, 1))
# Save as .png
ggsave(Figure1_combined,
filename = paste0(figurePath, "Fig_1_lower_part.png"),
dpi = 1000,
width = 170,
height = 90,
units = "mm")
# Save as .svg
ggsave(Figure1_combined,
filename = paste0(figurePath, "Fig_1_lower_part.svg"),
dpi = 1000,
width = 170,
height = 90,
units = "mm")
Descriptive statistics of exploration
Mean and SD of Raleigh vector length
# Report mean
str1 <- mean_SD_str2(sector_rayleigh$resultant_length, report_type, digits1, rounding_type)- Raleigh vector length: M = 0.37 (SD = 0.23)
Calculating the number of visits per participant
# Calculate the number of visited sectors
num_visitedSectors <- ddply(no_control, c("ppid"), summarise,
number = length_uniq(sector))
# Make report string
val <- num_visitedSectors$number
str1 <- mean_SD_str2(val, report_type, digits1, rounding_type)Number of visited sector per participant: M = 50 (SD = 2)
Describe two measures on which the novelty score is based
# Load image generated when creating the event files
load("event_tables/images/SpaNov_event_file_gradients.RData")
# Convert list to data frame
data <- as.data.frame(rbindlist(subj_list, idcol = "subject"))
# Add subjects' R number
data$ppid <- subjIDs_R[data$subject]
# Function to match runStartTime
find_runStartTime <- function(ppid, run){
# Get corresponding to find the run in the trial data
anonKey <- lookupTable$anonKey[lookupTable$Rnum == ppid[1]]
# Use the anonKey & run to get runStartTime
bool_index <- OLM_7T_trial_results$subject == anonKey &
OLM_7T_trial_results$block_num == run
runStartTime <- OLM_7T_trial_results$runStartTime[bool_index]
return(runStartTime[1])
}
# Use the function to find the correct run start time
data <- ddply(data, c("subject", "ppid", "run"), mutate,
runStartTime = find_runStartTime(ppid, run))
# Make time relative to the start of the run. The real onset times will be
# slightly different but this will not matter for this.
data$onset_rel <- data$onset - data$runStartTime
# Add run type
data$runType <- "encoding"
data$runType[data$run == 2 | data$run == 4] <- "retrieval"
# Subset to only encoding
data_sub <- data[data$runType == "encoding", ]
# Change run number to match the description in paper. Run 3 is Encoding Run 2
data_sub$run[data_sub$run == 3] <- 2
# Convert run to factor
data_sub$f_run <- as.factor(data_sub$run)
# Calculate descriptive stats for the measures
data_sub_agg1 <- ddply(data_sub, c("subject"), summarise,
max_visits = max(visits),
min_lastVisit = min(lastVisit, na.rm = TRUE),
max_lastVisit = max(lastVisit, na.rm = TRUE),
median_lastVisit = median(lastVisit, na.rm = TRUE),
mean_lastVisit = mean(lastVisit, na.rm = TRUE))
# Calculate the average duration of events from lastVisits for that exclude NA values
data_sub_agg2 <- ddply(na.omit(data_sub), c("subject"), summarise,
mean_duration = mean(duration),
median_duration = median(duration),
min_duration = min(duration),
max_duration = max(duration))
# Create strings for report
val <- data_sub_agg1$max_visits
str1 <- mean_SD_str2(val, report_type, digits1, rounding_type)
val <- data_sub_agg1$min_lastVisit
str2 <- mean_SD_str2(val, report_type, digits1, rounding_type, "s")
val <- data_sub_agg1$max_lastVisit
str3 <- mean_SD_str2(val, report_type, digits1, rounding_type, "s")
val <- data_sub_agg2$mean_duration
str4 <- mean_SD_str2(val, report_type, digits1, rounding_type, "s")
val <- data_sub_agg2$min_duration
str5 <- mean_SD_str2(val, report_type, digits1, rounding_type, "s")
val <- data_sub_agg2$max_duration
str6 <- mean_SD_str2(val, report_type, digits1, rounding_type, "s")- Maximum visit per sector for the participants was M = 18 (SD = 1.6) by the time of the second run.
- The minimum time between visits varied between participants M = 6.7 s (SD = 2.3 s).
- The maximum time between visits varied between participants M = 1300 s (SD = 310 s).
- The average duration of the events for the time since varied between participants M = 2.5 s (SD = 0.4 s) with an average minimum of M = 0.6 s (SD = 0.0023 s) and an average maximum of M = 12 s (SD = 5.3 s).
Locomotion during exploration
Report the average values for the three locomotion states
val <- locomotion_per$translation * 100
str1 <- mean_SD_str2(val, report_type, digits1, rounding_type, "%")
val <- locomotion_per$rotation * 100
str2 <- mean_SD_str2(val, report_type, digits1, rounding_type, "%")
val <- locomotion_per$stationary * 100
str3 <- mean_SD_str2(val, report_type, digits1, rounding_type, "%")- Translation: M = 55 % (SD = 8.2 %)
- Rotation: M = 24 % (SD = 4 %)
- Stationary: M = 21 % (SD = 6.3 %)
Travelled path length
# Create the report string
str1 <- mean_SD_str2(total_path_lengths$distance,
report_type, digits1, rounding_type, "vm")The average total path length was M = 3200 vm (SD = 140 vm)
Event duration
# Make report string
val <- event_df$duration
str1 <- mean_SD_str2(val, report_type, digits1, rounding_type)Event duration: M = 2.5 (SD = 2.1)
Bayesian hierarchical modelling of novelty score
# Get marginal results
mae_respone_noveltyScore_run1 <- avg_comparisons(m_noveltyScore_run1, variables = "s_onset_rel",
type = "response", re_formula = NULL)
mae_respone_noveltyScore_run2 <- avg_comparisons(m_noveltyScore_run2, variables = "s_onset_rel",
type = "response", re_formula = NULL)
str1 <- create_str_from_avg_comparisons(mae_respone_noveltyScore_run1, "SDs")
str2 <- create_str_from_avg_comparisons(mae_respone_noveltyScore_run2, "SDs")- Average marginal effect in Run 1: -0.48 SDs (95 % CI [-0.55, -0.41])
- Average marginal effect in Run 2: -0.71 SDs (95 % CI [-0.75, -0.66])
\[ \begin{align*} y & \sim \operatorname{Student}(\nu, \mu, \sigma) \\ \mu & = time + (time | subject) \end{align*} \]
Alternative way to describe: https://rpsychologist.com/r-guide-longitudinal-lme-lmer
Mapping spatial novelty across the hippocampal long axis
Load the results from the linear contrast (Level 1 > Level 2 > Level 3 > …). All xiis from this contrast start with GLM1.
- Positive values: Higher activity for familiar sectors (c1)
- Negative values: Higher activity for novel sectors (c2)
# Folder where the images are
ciftiFolder <- "/SpaNov/OLMe_7T_SpaNov_gradient_6lvl_cue-delay_smo4_MSMAll/cope7.feat/stats/vwc/"
# Other parameters
minClusterSize <- 20
# Choose the files
zMap_file <- paste0(path2imaging_results2, ciftiFolder, "results_lvl2cope1_dat_ztstat_c1.dscalar.nii")
pMap1_file <- paste0(path2imaging_results2, ciftiFolder, "results_lvl2cope1_dat_ztstat_cfdrp_c1.dscalar.nii")
pMap2_file <- paste0(path2imaging_results2, ciftiFolder, "results_lvl2cope1_dat_ztstat_cfdrp_c2.dscalar.nii")
clusterMap1_file <- paste0(path2imaging_results2, ciftiFolder, "results_lvl2cope1_dat_ztstat_cfdrp_c1_clusters.dscalar.nii")
clusterMap2_file <- paste0(path2imaging_results2, ciftiFolder, "results_lvl2cope1_dat_ztstat_cfdrp_c2_clusters.dscalar.nii")
betaMap_file <- paste0(path2imaging_results2, ciftiFolder, "Y1.dtseries.nii")
# Load all maps as xiis
GLM1_zMap_xii <- read_cifti(zMap_file, brainstructures = "all")
GLM1_pMap1_xii <- read_cifti(pMap1_file, brainstructures = "all")
GLM1_pMap2_xii <- read_cifti(pMap2_file, brainstructures = "all")
GLM1_clusterMap1_xii <- read_cifti(clusterMap1_file, brainstructures = "all")
GLM1_clusterMap2_xii <- read_cifti(clusterMap2_file, brainstructures = "all")
GLM1_betaMap_xii <- read_cifti(betaMap_file, brainstructures = "all")
# Create masks based on significant effects
HC_familiarity_mask <- GLM1_pMap1_xii$data$subcort > cutOff & str_detect(GLM1_pMap1_xii$meta$subcort$labels, pattern = "Hippocampus")
HC_novelty_mask <- GLM1_pMap2_xii$data$subcort > cutOff & str_detect(GLM1_pMap1_xii$meta$subcort$labels, pattern = "Hippocampus")
HC_familiarity_mask_L <- GLM1_pMap1_xii$data$subcort > cutOff & str_detect(GLM1_pMap1_xii$meta$subcort$labels, pattern = "Hippocampus-L")
HC_novelty_mask_L <- GLM1_pMap2_xii$data$subcort > cutOff & str_detect(GLM1_pMap1_xii$meta$subcort$labels, pattern = "Hippocampus-L")
HC_familiarity_mask_R <- GLM1_pMap1_xii$data$subcort > cutOff & str_detect(GLM1_pMap1_xii$meta$subcort$labels, pattern = "Hippocampus-R")
HC_novelty_mask_R <- GLM1_pMap2_xii$data$subcort > cutOff & str_detect(GLM1_pMap1_xii$meta$subcort$labels, pattern = "Hippocampus-R")
WB_familiarity_mask <- GLM1_pMap1_xii > cutOff
WB_novelty_mask <- GLM1_pMap2_xii > cutOff
# Extract & average values
HC_familiarity_values <- GLM1_betaMap_xii$data$subcort[HC_familiarity_mask, ]
HC_familiarity_values <- colMeans(HC_familiarity_values)
HC_novelty_values <- GLM1_betaMap_xii$data$subcort[HC_novelty_mask, ]
HC_novelty_values <- colMeans(HC_novelty_values)
HC_familiarity_values_L <- GLM1_betaMap_xii$data$subcort[HC_familiarity_mask_L, ]
HC_familiarity_values_L <- colMeans(HC_familiarity_values_L)
HC_familiarity_values_R <- GLM1_betaMap_xii$data$subcort[HC_familiarity_mask_R, ]
HC_familiarity_values_R <- colMeans(HC_familiarity_values_R)
HC_novelty_values_R <- GLM1_betaMap_xii$data$subcort[HC_novelty_mask_R, ]
HC_novelty_values_R <- colMeans(HC_novelty_values_R)
WB_novelty_values <- as.matrix(GLM1_betaMap_xii)
WB_novelty_values <- WB_novelty_values[as.matrix(WB_novelty_mask) == 1, ]
WB_novelty_values <- colMeans(WB_novelty_values)
WB_familiarity_values <- as.matrix(GLM1_betaMap_xii)
WB_familiarity_values <- WB_familiarity_values[as.matrix(WB_familiarity_mask) == 1, ]
WB_familiarity_values <- colMeans(WB_familiarity_values)
# Create cluster tables based on the maps
GLM1_cluster1 <- cifti_cluster_report(zMap_file,
clusterMap1_file,
surfLeft,
surfRight,
parcellationFile,
minClusterSize,
FALSE)GLM1_cluster2 <- cifti_cluster_report(zMap_file,
clusterMap2_file,
surfLeft,
surfRight,
parcellationFile,
minClusterSize,
FALSE)# Round values to in order to report them
GLM1_cluster1$cluster_values$zValue_max <- signif(GLM1_cluster1$cluster_values$zValue_max, digits1)
GLM1_cluster1$cluster_values$zValue_min <- signif(GLM1_cluster1$cluster_values$zValue_min, digits1)
GLM1_cluster1$cluster_values$cluster_mass <- signif(GLM1_cluster1$cluster_values$cluster_mass, digits1)
GLM1_cluster2$cluster_values$zValue_max <- signif(GLM1_cluster2$cluster_values$zValue_max, digits1)
GLM1_cluster2$cluster_values$zValue_min <- signif(GLM1_cluster2$cluster_values$zValue_min, digits1)
GLM1_cluster2$cluster_values$cluster_mass <- signif(GLM1_cluster2$cluster_values$cluster_mass, digits1)
# Write tsv files so it's easier to edit
## Positive clusters
### Combine to one table
region_labels <- GLM1_cluster1$cluster_labels[, 4]
combined_table <- cbind(GLM1_cluster1$cluster_values, region_labels)
### Add places between the region labels
combined_table$region_labels <- str_replace_all(combined_table$region_labels,
pattern = ",",
replacement = ", ")
### Write as .txt file
write.table(combined_table, file = "figures/SpaNov/tables/GLM1_c1.txt",
quote = FALSE,
row.names = FALSE,
sep = '\t')
## Negative clusters
### Combine to one table
region_labels <- GLM1_cluster2$cluster_labels[, 4]
combined_table <- cbind(GLM1_cluster2$cluster_values, region_labels)
### Add places between the region labels
combined_table$region_labels <- str_replace_all(combined_table$region_labels,
pattern = ",",
replacement = ", ")
### Write as .txt file
write.table(combined_table, file = "figures/SpaNov/tables/GLM1_c2.txt",
quote = FALSE,
row.names = FALSE,
sep = '\t')Create unthresholded hippocampal z-map
For unthresholded visualisation of the hippocampal GLM results, we create a version of the z-map that set everything other than the hippocampus to zero.
# (This chunk only needs to run once, so eval is set to FALSE)
# Extreme comparison: Level 1 vs. Level 6
ciftiFolder <- "/SpaNov/OLMe_7T_SpaNov_gradient_6lvl_cue-delay_smo4_MSMAll/cope14.feat/stats/vwc/"
zMap_file1 <- paste0(path2imaging_results2, ciftiFolder, "results_lvl2cope1_dat_ztstat_c1.dscalar.nii")
zMap_xii1 <- read_cifti(zMap_file1, brainstructures = "all",
surfL_fname = surfLeft, surfR_fname = surfRight)
# Linear contrast from Level 1 to Level 6
ciftiFolder <- "/SpaNov/OLMe_7T_SpaNov_gradient_6lvl_cue-delay_smo4_MSMAll/cope7.feat/stats/vwc/"
zMap_file2 <- paste0(path2imaging_results2, ciftiFolder, "results_lvl2cope1_dat_ztstat_c1.dscalar.nii")
zMap_xii2 <- read_cifti(zMap_file2, brainstructures = "all",
surfL_fname = surfLeft, surfR_fname = surfRight)
# Create mask for everything other than right and left hippocampus
currentMask <- zMap_xii1$meta$subcort$labels != "Hippocampus-R" &
zMap_xii1$meta$subcort$labels != "Hippocampus-L"
# Create new variables
zMap_xii1_HC <- zMap_xii1
zMap_xii2_HC <- zMap_xii2
# Set everything within this mask to zero
zMap_xii1_HC$data$subcort[currentMask, ] <- 0
zMap_xii2_HC$data$subcort[currentMask, ] <- 0
# Also set cortical values to zero
zMap_xii1_HC$data$cortex_left <- matrix(as.integer(0), nrow = 29696, ncol = 1)
zMap_xii1_HC$data$cortex_right <- matrix(as.integer(0), nrow = 29716, ncol = 1)
zMap_xii2_HC$data$cortex_left <- matrix(as.integer(0), nrow = 29696, ncol = 1)
zMap_xii2_HC$data$cortex_right <- matrix(as.integer(0), nrow = 29716, ncol = 1)
# Create new names
new_zMap_file1 <- str_replace(zMap_file1, pattern = ".dscalar.nii",
replacement = "_onlyHC.dscalar.nii")
new_zMap_file2 <- str_replace(zMap_file2, pattern = ".dscalar.nii",
replacement = "_onlyHC.dscalar.nii")
# Write new cifti files
write_cifti(xifti = zMap_xii1_HC, cifti_fname = new_zMap_file1)
write_cifti(xifti = zMap_xii2_HC, cifti_fname = new_zMap_file2)Spatial novelty-gradient in hippocampus
Calculate GS & tSNR in hippocampus
# Load tSNR and GS values
load("data/ignore_fMRI_version1/tSNR_GS_maps.RData")
# Add file ID
GS_tSNR_df$rowID <- 1:nrow(GS_tSNR_df)
# Subset to subjects included in the analysis
GS_tSNR_df <- GS_tSNR_df[GS_tSNR_df$GS_subject %in% subjIDs_R, ]
# Calculate average GS & tSNR maps
GS_merged <- merge_xifti(xifti_list = GS_list[GS_tSNR_df$rowID])
GS_avg <- apply_xifti(GS_merged, margin = 1, mean)
tSNR_merged <- merge_xifti(xifti_list = tSNR_list[GS_tSNR_df$rowID])
tSNR_avg <- apply_xifti(tSNR_merged, margin = 1, mean)
# Get MNI coordinates
HC_data_L <- MNI_coord[str_starts(MNI_coord$region, pattern = "Hippocampus-L"), ]
HC_data_R <- MNI_coord[str_starts(MNI_coord$region, pattern = "Hippocampus-R"), ]
HC_data <- rbind(HC_data_L, HC_data_R)
# Get bool index
HC_index_L <- GS_list[[1]]$meta$subcort$labels == "Hippocampus-L"
HC_index_R <- GS_list[[1]]$meta$subcort$labels == "Hippocampus-R"
# Results list
GS_HC_list <- list()
tSNR_HC_list <- list()
# Extract hippocampus values by loop over all subjects
for(i in 1:length(subjIDs_R)){
# GS
## Extract values from Run 1 & 2
run1_xii <- GS_list[[GS_tSNR_df$rowID[GS_tSNR_df$GS_run == "RUN1" & GS_tSNR_df$GS_subject == subjIDs_R[i]]]]
run2_xii <- GS_list[[GS_tSNR_df$rowID[GS_tSNR_df$GS_run == "RUN2" & GS_tSNR_df$GS_subject == subjIDs_R[i]]]]
run1_HC <- c(run1_xii$data$subcort[HC_index_L], run1_xii$data$subcort[HC_index_R])
run2_HC <- c(run2_xii$data$subcort[HC_index_L], run2_xii$data$subcort[HC_index_R])
## Average across runs
avg_HC <- (run1_HC + run2_HC)/2
## Add to data.frame
temp_df <- HC_data
temp_df$subject <- subjIDs_R[i]
temp_df$GS <- avg_HC
GS_HC_list[[i]] <- temp_df
# GS
## Extract values from Run 1 & 2
run1_xii <- tSNR_list[[GS_tSNR_df$rowID[GS_tSNR_df$GS_run == "RUN1" & GS_tSNR_df$GS_subject == subjIDs_R[i]]]]
run2_xii <- tSNR_list[[GS_tSNR_df$rowID[GS_tSNR_df$GS_run == "RUN2" & GS_tSNR_df$GS_subject == subjIDs_R[i]]]]
run1_HC <- c(run1_xii$data$subcort[HC_index_L], run1_xii$data$subcort[HC_index_R])
run2_HC <- c(run2_xii$data$subcort[HC_index_L], run2_xii$data$subcort[HC_index_R])
## Average across runs
avg_HC <- (run1_HC + run2_HC)/2
## Add to data.frame
temp_df <- HC_data
temp_df$subject <- subjIDs_R[i]
temp_df$tSNR <- avg_HC
tSNR_HC_list[[i]] <- temp_df
}
# Convert to data frame
GS_HC_df <- rbindlist(GS_HC_list)
tSNR_HC_df <- rbindlist(tSNR_HC_list)
GS_HC_df <- as.data.frame(GS_HC_df)
tSNR_HC_df <- as.data.frame(tSNR_HC_df)
# Simple A vs. P comparison
## Divide into A & P
GS_HC_df$AP <- ifelse(GS_HC_df$y > -21, "Anterior", "Posterior")
tSNR_HC_df$AP <- ifelse(tSNR_HC_df$y > -21, "Anterior", "Posterior")
## Calculate average for each participant
GS_HC_df_AP <- ddply(GS_HC_df, c("subject", "AP"), summarise, GS = mean(GS))
tSNR_HC_df_AP <- ddply(tSNR_HC_df, c("subject", "AP"), summarise, tSNR = mean(tSNR))
## Run tests
### GS
GS_A <- GS_HC_df_AP$GS[GS_HC_df_AP$AP == "Anterior"]
GS_P <- GS_HC_df_AP$GS[GS_HC_df_AP$AP == "Posterior"]
diff <- GS_A - GS_P
GS_d <- round(mean(diff)/sd(diff), 2)
GS_p <- t.test(diff)$p.value
GS_BF10 <- signif(reportBF(ttestBF(diff)), 3)
### tSNR
tSNR_A <- tSNR_HC_df_AP$tSNR[tSNR_HC_df_AP$AP == "Anterior"]
tSNR_P <- tSNR_HC_df_AP$tSNR[tSNR_HC_df_AP$AP == "Posterior"]
diff <- tSNR_A - tSNR_P
tSNR_d <- round(mean(diff)/sd(diff), 2)
tSNR_p <- t.test(diff)$p.value
tSNR_BF10 <- signif(reportBF(ttestBF(diff)), 3)
# Calculate percentile of averge hippocampus tSNR
tSNR_A_percentile <- round(mean(as.matrix(tSNR_avg) < mean(tSNR_A)), 2)
tSNR_P_percentile <- round(mean(as.matrix(tSNR_avg) < mean(tSNR_P)), 2)
# Get values from the XIFTIs
HC_tSNR_avg <- c(tSNR_avg$data$subcort[HC_index_L], tSNR_avg$data$subcort[HC_index_R])
HC_GS_avg <- c(GS_avg$data$subcort[HC_index_L], GS_avg$data$subcort[HC_index_R])
## Create a data frame
HC_tSNR_GS_df <- data.frame(tSNR = HC_tSNR_avg, GS = HC_GS_avg,
region = rep(c("Hippocampus-L", "Hippocampus-R"),
times = c(sum(HC_index_L), sum(HC_index_R))))Permutation analyses
# Load hippocampal gradient data
load("intermediate_data/SpaNov_gradient_data_cue-delay.RData")
# Add tSNR & GS
HC_data$tSNR <- HC_tSNR_GS_df$tSNR
HC_data$GS <- HC_tSNR_GS_df$GS
# In 5b1e712, carefully check that the order HC_data and HC_tSNR_GS_df is
# the same. Specifically I calculated the position values for HC_tSNR_GS_df
# and then checked if mean(HC_data$position == HC_tSNR_GS_df$position) == 1
# Average across the AP position
HC_data_agg_pos <- ddply(HC_data, c("position", "Hemisphere"), summarise,
n = length(min), min = mean(min), max = mean(max),
tSNR = mean(tSNR), GS = mean(GS))
# Average across position and hemisphere
HC_data_agg_pos2 <- ddply(HC_data, c("position"), summarise,
n = length(min), min = mean(min), max = mean(max),
tSNR = mean(tSNR), GS = mean(GS))
# Labels
conN <- 6
novLabels <- paste0('Level ', 1:conN)permutation_analysis <- function(data, lm_formula, nIter, colName, imageName){
# Initialise results list
results <- list()
# Select correct column for analysis and create new data frame
data$val <- data[, colName]
data2shuffle <- data
# Calculate empirical values and save to list
results$lm <- lm(lm_formula, data = data)
numCoef <- length(results$lm$coefficients) - 1 # Ignoring the intercept
# Start cluster
my.cluster <- parallel::makeCluster(detectCores() - 2, type = "PSOCK")
#register it to be used by %dopar%
doParallel::registerDoParallel(cl = my.cluster)
# Run parallel loop
permuted_values <- foreach(i = 1:nIter, .combine = 'c', .packages = 'plyr') %dopar% {
data2shuffle$val <- sample(data$val)
# Fit model
temp_lm <- lm(lm_formula, data = data2shuffle)
# Add values
temp_est <- as.data.frame(matrix(as.numeric(temp_lm$coefficients)[-1], ncol = numCoef))
names(temp_est) <- names(results$lm$coefficients)[-1]
list(temp_est)
}
# Stop cluster again
parallel::stopCluster(cl = my.cluster)
# Add to results
results$permuted_values <- as.data.frame(rbindlist(permuted_values))
# Save to disk
if(!missing(imageName)){
save(results, file = paste0("intermediate_data/", imageName))
}
# Return value
return(results)
}Minimum
# Rename to be unique
grad_HC1 <- HC_data_agg_pos
# Check if the code needs to be run again. This is done by checking the md5 has
# for the data frame that is used in the calculation if it is different from the
# hash sum saved in md5_hash_table.csv, it is re-run. This is to avoid having to
# re-run everything each time I work on the data.
runCodeAgain1 <- check_if_md5_hash_changed(grad_HC1, hash_table_name = "SpaNov_md5_hash_table.csv")Model with interaction with hemisphere
# Seed
set.seed(19911225)
# Other input
lm_formula <- "val ~ position * Hemisphere + tSNR + GS"
nIter <- 100000
colName <- "min"
imageName <- "SpaNov_permut_HC_grad_analysis1_cue-delay.RData"
# Run if necessary
if(runCodeAgain1){
grad_min_permut1 <- permutation_analysis(grad_HC1, lm_formula,
nIter, colName, imageName)
} else {
load(paste0("intermediate_data/", imageName))
grad_min_permut1 <- results
}
# Select values for plotting
dist <- grad_min_permut1$permuted_values[,5]
critVal <- grad_min_permut1$lm$coefficients[6]
grad_min_permut1_p5 <- pValue_from_nullDist(critVal, dist, "two.sided")Left hippocampus
# Seed
set.seed(19911225)
# Other input
lm_formula <- "val ~ position + tSNR + GS"
nIter <- 100000
colName <- "min"
imageName <- "SpaNov_permut_HC_grad_analysis1_cue-delay_L.RData"
# Run if necessary
if(runCodeAgain1){
grad_min_permut1_L <- permutation_analysis(grad_HC1[grad_HC1$Hemisphere == 'left', ], lm_formula,
nIter, colName, imageName)
} else {
load(paste0("intermediate_data/", imageName))
grad_min_permut1_L <- results
}
# Select values for plotting
dist <- grad_min_permut1_L$permuted_values[,1]
critVal <- grad_min_permut1_L$lm$coefficients[2]
grad_min_permut1_L_p1 <- pValue_from_nullDist(critVal, dist, "two.sided")Right hippocampus
# Seed
set.seed(19911225)
# Other input
lm_formula <- "val ~ position + tSNR + GS"
nIter <- 100000
colName <- "min"
imageName <- "SpaNov_permut_HC_grad_analysis1_cue-delay_R.RData"
# Run if necessary
if(runCodeAgain1){
grad_min_permut1_R <- permutation_analysis(grad_HC1[grad_HC1$Hemisphere == 'right', ], lm_formula,
nIter, colName, imageName)
} else {
load(paste0("intermediate_data/", imageName))
grad_min_permut1_R <- results
}
# Select values for plotting
dist <- grad_min_permut1_R$permuted_values[,1]
critVal <- grad_min_permut1_R$lm$coefficients[2]
grad_min_permut1_R_p1 <- pValue_from_nullDist(critVal, dist, "two.sided")Maximum
Model with interaction with hemisphere
# Seed
set.seed(19911225)
# Other input
lm_formula <- "val ~ position * Hemisphere + tSNR + GS"
nIter <- 100000
colName <- "max"
imageName <- "SpaNov_permut_HC_grad_analysis3_cue-delay.RData"
# Run if necessary
if(runCodeAgain1){
grad_max_permut1 <- permutation_analysis(grad_HC1, lm_formula,
nIter, colName, imageName)
} else {
load(paste0("intermediate_data/", imageName))
grad_max_permut1 <- results
}
# Select values for plotting
dist <- grad_max_permut1$permuted_values[,5]
critVal <- grad_max_permut1$lm$coefficients[6]
grad_max_permut1_p5 <- pValue_from_nullDist(critVal, dist, "two.sided")Average across hemispheres
# Rename to be unique
grad_HC2 <- HC_data_agg_pos2
# Check if the code needs to be run again. This is done by checking the md5 has
# for the data frame that is used in the calculation if it is different from the
# hash sum saved in md5_hash_table.csv, it is re-run. This is to avoid having to
# re-run everything each time I work on the data.
runCodeAgain2 <- check_if_md5_hash_changed(grad_HC2, hash_table_name = "SpaNov_md5_hash_table.csv")
# Seed
set.seed(20131)
# Other input
lm_formula <- "val ~ position + tSNR + GS"
nIter <- 100000
colName <- "max"
imageName <- "SpaNov_permut_HC_grad_analysis4_cue-delay.RData"
# Run if necessary
if(runCodeAgain2){
grad_max_permut1_avg <- permutation_analysis(grad_HC2, lm_formula,
nIter, colName, imageName)
} else {
load(paste0("intermediate_data/", imageName))
grad_max_permut1_avg <- results
}
# Select values for plotting
dist <- grad_max_permut1_avg$permuted_values[,1]
critVal <- grad_max_permut1_avg$lm$coefficients[2]
grad_max_permut1_avg_p1 <- pValue_from_nullDist(critVal, dist, "two.sided")Report stat of permutation analysis
# Calculate effect sizes
grad_min_permut1_eff <- eta_squared(car::Anova(grad_min_permut1$lm, type = 2))
grad_min_permut1_L_eff <- eta_squared(car::Anova(grad_min_permut1_L$lm, type = 2))
grad_min_permut1_R_eff <- eta_squared(car::Anova(grad_min_permut1_R$lm, type = 2))
grad_max_permut1_R_eff <- eta_squared(car::Anova(grad_max_permut1$lm, type = 2))
grad_max_permut1_avg_eff <- eta_squared(car::Anova(grad_max_permut1_avg$lm, type = 2))
# Create data frame
## Create variables
tmp_formula <- c("min ~ position * Hemisphere + tSNR + GS", "min (left) ~ position + tSNR + GS",
"min (right) ~ position + tSNR + GS", "max ~ position * Hemisphere + tSNR + GS",
"max (avg) ~ position + tSNR + GS")
tmp_coef_name <- c("interaction", "position", "position", "interaction", "position")
tmp_coef_val <- c(grad_min_permut1$lm$coefficients[6], grad_min_permut1_L$lm$coefficients[2],
grad_min_permut1_R$lm$coefficients[2], grad_max_permut1$lm$coefficients[6],
grad_max_permut1_avg$lm$coefficients[2])
tmp_coef_es <- c(grad_min_permut1_eff$Eta2_partial[5], grad_min_permut1_L_eff$Eta2_partial[1],
grad_min_permut1_R_eff$Eta2_partial[1], grad_max_permut1_R_eff$Eta2_partial[5],
grad_max_permut1_avg_eff$Eta2_partial[1])
tmp_coef_p <- c(grad_min_permut1_p5, grad_min_permut1_L_p1, grad_min_permut1_R_p1,
grad_max_permut1_p5, grad_max_permut1_avg_p1)
## Convert to numeric values and flip
tmp_coef_val <- -as.numeric(tmp_coef_val)
## Add to one data frame
tab_df <- data.frame(Formula = tmp_formula, Coefficient = tmp_coef_name,
beta = tmp_coef_val, eta_squared = tmp_coef_es, p = tmp_coef_p)
## Round columns
tab_df$beta <- round(tab_df$beta, 3)
tab_df$eta_squared <- round(tab_df$eta_squared, 3)
## Create p-values by looping over all values
for(i in seq_along(tab_df$p)){
tab_df$p[i] <-paste("p", pValue(as.numeric(tab_df$p[i])))
}
# Show table
kable(tab_df)| Formula | Coefficient | beta | eta_squared | p |
|---|---|---|---|---|
| min ~ position * Hemisphere + tSNR + GS | interaction | 0.022 | 0.030 | p = .021 |
| min (left) ~ position + tSNR + GS | position | 0.021 | 0.053 | p = .021 |
| min (right) ~ position + tSNR + GS | position | 0.039 | 0.154 | p < .001 |
| max ~ position * Hemisphere + tSNR + GS | interaction | 0.013 | 0.011 | p = .194 |
| max (avg) ~ position + tSNR + GS | position | 0.031 | 0.104 | p < .001 |
Indidiviudal slopes
# Load individual slopes
load("data/ignore_fMRI_version1/extracted_values/SpavNov_gradient_subject-slopes_20250721_073753.RData")
# Add significance
subjSlopes$sig <- subjSlopes$`Pr(>|t|)` < .05
# Individual slopes from gradient analysis
gradient_subject_slopes_hemisphere <- -subjSlopes$Estimate[subjSlopes$type == "position" & subjSlopes$effect == "position:Hemisphereright"]
gradient_subject_slopes_position_L <- -subjSlopes$Estimate[subjSlopes$type == "position_L" & subjSlopes$effect == "position"]
gradient_subject_slopes_position_R <- -subjSlopes$Estimate[subjSlopes$type == "position_R" & subjSlopes$effect == "position"]
gradient_subject_slopes_position_L_sig <- subjSlopes$sig[subjSlopes$type == "position_L" & subjSlopes$effect == "position"]
gradient_subject_slopes_position_R_sig <- subjSlopes$sig[subjSlopes$type == "position_R" & subjSlopes$effect == "position"]
# Effect of hemispheres
x <- gradient_subject_slopes_hemisphere
d1 <- signif(mean(x)/sd(x), digits1)
str1 <- mean_SD_str2(x, report_type, digits1, rounding_type)
BF10_1 <- reportBF(ttestBF(x))
# Effect in left hemisphere
x <- gradient_subject_slopes_position_L
d3 <- signif(mean(x)/sd(x), digits1)
str3 <- mean_SD_str2(x, report_type, digits1, rounding_type)
BF10_3 <- reportBF(ttestBF(x))
# Effect in right hemisphere
x <- gradient_subject_slopes_position_R
d4 <- signif(mean(x)/sd(x), digits1)
str4 <- mean_SD_str2(x, report_type, digits1, rounding_type)
BF10_4 <- reportBF(ttestBF(x))In the right hippocampus, individual slopes were different from zero,M = 0.0093 (SD = 0.017), BF10 = 148.74, d = 0.54, while they was weak evidence for that null hypothesis in the left hippocampus,M = 0.0037 (SD = 0.017), BF10 = 0.49, d = 0.21. However, evidence for a difference between hemisphere in terms of individual slopes was inconclusive,M = 0.0056 (SD = 0.017), BF10 = 2.04, d = 0.32.
Figure 2
# Create new data frame just for plotting
HC_data_plotting <- HC_data_agg_pos
HC_data_plotting$Hemisphere2 <- ifelse(HC_data_plotting$Hemisphere == "right",
"Right hippocampus", "Left hippocampus")
# Write ,csv file
write.csv(x = HC_data_plotting, file = "other_stuff/source_data_files/HC_scatter_plot.csv",
quote = FALSE, row.names = FALSE)
# Create Figure 2
## Create minimum scatter plot for both hemispheres
scatter_min <- ggplot(HC_data_plotting, aes(x = -position, y = min)) +
facet_grid(~Hemisphere2, scales = "free_x") +
geom_point(aes(colour = min)) +
geom_smooth(method = "lm", formula = y ~ x, colour = "black") +
labs(x = "Distance to most anterior part (mm)", y = "Novelty\npreference") +
scale_x_continuous(breaks = seq(from = -45, to = 0, by = 15), labels = c("45", "30", "15", "0")) +
scale_y_continuous(breaks = 1:6, labels = paste("Lvl", 1:6), limits = c(1, 6)) +
scale_colour_viridis_c(option = "H", limits = c(1, 6), direction = -1) +
theme_journal() +
theme(strip.background = element_rect(color="white", fill="white"),
axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
panel.spacing.x = unit(7.5, "mm")) +
theme(legend.position = "none")
## Create plots for individual slopes
### Create data frame for plotting
individual_slope_df <- data.frame(Effect = rep(c("Right hippocampus", "Left hippocampus", "Interaction"), each = 56),
Coefficient = c(gradient_subject_slopes_position_R,
gradient_subject_slopes_position_L,
gradient_subject_slopes_hemisphere))
# Write ,csv file
write.csv(x = individual_slope_df, file = "other_stuff/source_data_files/HC_individual_slopes.csv",
quote = FALSE, row.names = FALSE)
### Create plot
individual_slopes_plot <- ggplot(individual_slope_df, aes(x = Effect, y = Coefficient)) +
geom_hline(yintercept = 0, linetype = 2) +
geom_jitter(colour = boxplot_point_col, width = 0.25) +
geom_boxplot(colour = boxplot_border_col, alpha = 0.7, outlier.shape = NA, width = 0.4,
fill = base_col[2]) +
stat_summary(geom = "point", fun = "mean", col = 'white', shape = 24,
fill = base_col[2], position = position_dodge(width = 0.75), size = 0.8) +
theme_journal() +
theme(axis.title.y = element_text(margin = margin(t = 0, r = 4.5, b = 0, l = 0))) +
labs(title = "", x = "", y = "Individual\nregression coefficients")
## Combine both plots
p <- plot_grid(scatter_min, individual_slopes_plot, ncol = 1, align = "hv", axis = "tblr")
## Save as .png
ggsave(p, filename = paste0(figurePath, "Fig_2c_scattter_and_boxplot.png"),
dpi = dpi, width = 99, height = 83.5, units = "mm")
## Save as .svg
ggsave(p, filename = paste0(figurePath, "Fig_2c_scattter_and_boxplot.svg"),
dpi = dpi, width = 99, height = 83.5, units = "mm")
Extension of spatial novelty gradient to the posterior medial cortex
Figure 3
Yeo 7 parcellation ridge plot
# Get the names of the networks (-1 to remove the first row that just contains ???)
Yeo7_labels <- Yeo7_xii$meta$cifti$labels$parcels[-1, ]
Yeo7_labels$Region <- row.names(Yeo7_labels)
row.names(Yeo7_labels) <- NULL
# Convert the RGB values to hexcodes
Yeo7_labels$hexcol <- rgb(Yeo7_labels$Red, Yeo7_labels$Green, Yeo7_labels$Blue,
maxColorValue = 1)
# Create basic data frame
Yeo_zValues_df <- data.frame(zValue = c(GLM1_zMap_xii$data$cortex_left, GLM1_zMap_xii$data$cortex_right),
Yeo7 = c(Yeo7_xii$data$cortex_left, Yeo7_xii$data$cortex_right))
# Add network label
Yeo_zValues_df$net_label <- NA
for(i in 1:nrow(Yeo7_labels)){
# Select all rows in Yeo_zValues_df that have this key
bool_index <- Yeo_zValues_df$Yeo7 == Yeo7_labels$Key[i]
Yeo_zValues_df$net_label[bool_index] <- Yeo7_labels$Region[i]
}
# Find lowest z value that's significant
## Get the p-values and the z-values
GLM1_pValues1 <- get_all_points_from_xifti(GLM1_pMap1_xii)
GLM1_pValues2 <- get_all_points_from_xifti(GLM1_pMap2_xii)
GLM1_zValues <- get_all_points_from_xifti(GLM1_zMap_xii)
## Subset to only significant values
GLM1_zValues1 <- GLM1_zValues[GLM1_pValues1 > cutOff]
GLM1_zValues2 <- GLM1_zValues[GLM1_pValues2 > cutOff]
# Exclude vertices that are not included in Yeo 7
Yeo_zValues_df_sub <- Yeo_zValues_df[Yeo_zValues_df$Yeo7 != 0, ]
# Order according to the mean of the network
Yeo_zValues_df_agg <- ddply(Yeo_zValues_df_sub, c("net_label"), summarise,
avg_z = mean(zValue))
## First order alphabetically so we can apply the same order to GLM1_Yeo7
Yeo_zValues_df_agg <- Yeo_zValues_df_agg[order(as.character(Yeo_zValues_df_agg$net_label)), ]
## Now order based on the mean
mean_order <- order(Yeo_zValues_df_agg$avg_z)
Yeo_zValues_df_agg <- Yeo_zValues_df_agg[mean_order, ]
# Order the colours in the same way
Yeo7_labels <- Yeo7_labels[order(as.character(Yeo7_labels$Region)), ]
Yeo7_labels <- Yeo7_labels[mean_order, ]
# Add the full names
Yeo7_labels$fullNames <- find_values_thru_matching(Yeo7_abbr, Yeo7_fullNames, Yeo7_labels$Region)
# Make net_label a factor with the correct order
Yeo_zValues_df_sub$net_label <- factor(Yeo_zValues_df_sub$net_label,
levels = Yeo_zValues_df_agg$net_label,
labels = Yeo7_labels$fullNames,
ordered = TRUE)
# Create plot
p1 <- ggplot(Yeo_zValues_df_sub, aes(x = zValue, y = net_label, fill = net_label)) +
geom_density_ridges(linewidth = 0.3) +
scale_fill_manual(values = Yeo7_labels$hexcol) +
geom_vline(xintercept = max(GLM1_zValues2), linetype = 2, linewidth = 0.3) +
geom_vline(xintercept = min(GLM1_zValues1), linetype = 2, linewidth = 0.3) +
theme_journal() +
coord_cartesian(ylim = c(0.5, 8.5)) +
theme(legend.position = "none", plot.margin = unit(c(1,1,1,1), "mm")) +
labs(x = "z-statistic", y = "")
# Save as .png
ggsave(p1,
filename = paste0(figurePath, "Fig_3a_Yeo7_ridge.png"),
dpi = dpi,
width = 60,
height = 30,
units = "mm")
# Save as .svg
ggsave(p1,
filename = paste0(figurePath, "Fig_3a_Yeo7_ridge.svg"),
dpi = dpi,
width = 60,
height = 30,
units = "mm")Margulies’ connectivity gradients
# Load all 10 gradients from brainstat
# Source: https://brainstat.readthedocs.io/en/master/index.html
nGrad <- 10
prefix <- "data/ignore_fMRI_version1/brainstat/Gradient_"
# Get the medial walls from donour
mwm_L <- Yeo7_xii$meta$cortex$medial_wall_mask$left
mwm_R <- Yeo7_xii$meta$cortex$medial_wall_mask$right
# Create data frame with
Margulies_grad <- as.data.frame(matrix(NA, nrow = 59412, ncol = 10))
names(Margulies_grad) <- paste0("Grad_", 1:10)
# Loop through all gradients
for(i in 1:nGrad){
# Create a new xifti by reading the gifti files
tmp_xii <- read_xifti2(cortexL = paste0(prefix, i, "_L.shape.gii"),
surfL = surfLeft,
mwall_values = c(0),
cortexR = paste0(prefix, i, "_R.shape.gii"),
surfR = surfRight)
# Use the medial wall values from donour and apply them to the loaded file.
tmp_xii$data$cortex_left[!mwm_L] <- NA
tmp_xii$data$cortex_right[!mwm_R] <- NA
tmp_xii <- move_to_mwall(tmp_xii, values = c(NA))
# Add to prepared data frame
Margulies_grad[, i] <- c(tmp_xii$data$cortex_left, tmp_xii$data$cortex_right)
}
# Add Yeo 7 to the data frame
Margulies_grad$Yeo7_name <- Yeo_zValues_df$net_label
# Remove the missing values
Margulies_grad_sub <- Margulies_grad[!is.na(Margulies_grad$Yeo7_name), ]
# Plot
p1 <- ggplot(Margulies_grad_sub, aes(x = Grad_2, y = Grad_1, colour = Yeo7_name)) +
geom_point(alpha = 0.2, size = 0.3, shape = 16) +
theme_void() +
theme(legend.position = "none") +
scale_colour_manual(values = Yeo7_colours)
# Save file
ggsave(p1,
filename = paste0(figurePath, "Fig_3a_Yeo7_Margulies_space.png"),
dpi = dpi,
width = 30,
height = 30,
units = "mm")
# Significant vertices for GLM1
Margulies_grad$GLM_1_direction <- NA
Margulies_grad$GLM_1_direction[c(GLM1_pMap1_xii$data$cortex_left, GLM1_pMap1_xii$data$cortex_right) > cutOff] <- "Familiarity"
Margulies_grad$GLM_1_direction[c(GLM1_pMap2_xii$data$cortex_left, GLM1_pMap2_xii$data$cortex_right) > cutOff] <- "Novelty"
# Create subset of only significant vertices
Margulies_grad_sub <- Margulies_grad[!is.na(Margulies_grad$GLM_1_direction), ]
# Plot
p1 <- ggplot(Margulies_grad_sub, aes(x = Grad_2, y = Grad_1, colour = GLM_1_direction)) +
geom_point(alpha = 0.2, size = 0.3, shape = 16) +
theme_void() +
scale_colour_manual(values = cool_warm_colours) +
theme(legend.position = "none")
# Save file
ggsave(p1,
filename = paste0(figurePath, "Fig_3a_NovFam_Margulies_space.png"),
dpi = dpi,
width = 70,
height = 70,
units = "mm")# Splitting the dataset into the Training set and Test set
# https://www.geeksforgeeks.org/classifying-data-using-support-vector-machinessvms-in-r/
svm_df <- data.frame(Grad_1 = Margulies_grad_sub$Grad_1,
Grad_2 = Margulies_grad_sub$Grad_2,
Modulation = ifelse(Margulies_grad_sub$GLM_1_direction == "Novelty", 1, 0))
# Seed for this analysis
set.seed(123)
# Split the data into training and test
split <- sample.split(svm_df$Modulation, SplitRatio = 0.75)
training_set <- subset(svm_df, split == TRUE)
test_set <- subset(svm_df, split == FALSE)
# Feature scaling
training_set[-3] <- scale(training_set[-3])
test_set[-3] <- scale(test_set[-3])
# Fitting SVM to the Training set
classifier <- svm(formula = Modulation ~ .,
data = training_set,
type = 'C-classification',
kernel = 'radial')
# Remove NA
test_set <- na.omit(test_set)
# Predicting the Test set results
y_pred <- predict(classifier, newdata = test_set)
cm <- table(test_set[, 3], y_pred)/ sum(table(test_set[, 3], y_pred))
cm <- round(cm *100, 2)
accuracy <- mean(test_set[, 3] == y_pred)The SVM-classification accuracy is 95.17% for significant novelty vs. familiarity vertices in Margulies’ et al. (2015) gradient space.
Resting-state connectivity between spatial novelty gradients in hippocampus and the posterior parietal cortex
Figure 4
############ connectivity_schematic
# Create hippocampus gradient bar
SC_gradient_xii <- read_cifti("cifti_results/SpaNov_gradient_wholebrain_min_cue-delay.dlabel.nii", brainstructures = c("subcort"))
# Create hippocampus only gradient map and write as CIFTI
HC_gradient_xii <- SC_gradient_xii
HC_gradient_xii$data$subcort[!str_detect(SC_gradient_xii$meta$subcort$labels, pattern = "Hippocampus"), 1] <- as.integer(0)
#write_cifti(HC_gradient_xii, cifti_fname = "cifti_results/SpaNov_gradient_HC_min_cue-delay.dlabel.nii")
############# Model results
# Add Marginal results
load("fitted_brms_models/SpaNov_m_conn2.Rdata")
# Rename model & set variable of interest
m <- m_conn2
voi <- "diagonal"
# Plot means of the conditions
## Use prediction
pred_means <- predictions(m, newdata = datagrid(diagonal = unique, f_gradient_level_cortex = unique,
f_gradient_level_HC = unique, subject = unique),
by = voi, re_formula = NULL, type = "response")
pred_means <- as.data.frame(pred_means)
## Get posterior draws
pred_means_draws <- get_draws(pred_means)
# Plot means for each participant
pred_means_subject <- avg_comparisons(m, variables = voi, by = c("subject"))
pred_means_subject <- as.data.frame(pred_means_subject)
## Get posterior draws
pred_means_subject_draws <- get_draws(pred_means_subject)
pred_means_subject_draws <- as.data.frame(pred_means_subject_draws)
# Order by mean difference
subject_ordered <- pred_means_subject$subject[order(pred_means_subject$estimate, decreasing = TRUE)]
pred_means_subject_draws$f_subject <- factor(pred_means_subject_draws$subject,
levels = subject_ordered, ordered = TRUE)
# Plot
# Width: a vector of probabilities to use that determine the widths of the resulting intervals
# Here that to 95%
pred_means_subject_draws_p <- ggplot(pred_means_subject_draws, aes(y = f_subject, x = -draw, colour = f_subject)) +
geom_vline(xintercept = 0, linetype = "dashed") +
stat_pointinterval(.width = 0.95, linewidth = 0.25, point_size = 0.05) +
labs(x = "Within novelty - between\nnovelty connectivity", y = "Individuals") +
theme_journal() +
theme(plot.title = element_text(hjust = 0.5),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
plot.subtitle = element_text(hjust = 0.5),
legend.position = "none")
# Average difference between the conditions
## Calculate marginal average effect is the contrast 0 vs. 1
mae1_respone <- avg_comparisons(m, variables = voi, type = "response", re_formula = NULL)
## Get posterior draws
mae_response_draws <- get_draws(mae1_respone)
## Flipping sign because it makes more sense
mae_response_draws$draw <- -mae_response_draws$draw
## Create plot
p1_mae_draws <- ggplot(mae_response_draws, aes(x = draw)) +
geom_vline(xintercept = 0, linetype = "dashed") +
stat_halfeye(fill = base_col[2], point_size = 1.5) +
labs(x = "Within novelty - between\nnovelty connectivity", y = "Density") +
theme_journal() +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
# Save as .png
ggsave(plot_grid(p1_mae_draws, pred_means_subject_draws_p, align = "h"),
filename = paste0(figurePath, "Fig_4_FRSC.png"),
dpi = dpi,
width = 90,
height = 45,
units = "mm")
# Save as .svg
ggsave(plot_grid(p1_mae_draws, pred_means_subject_draws_p, align = "h"),
filename = paste0(figurePath, "Fig_4_FRSC.svg"),
dpi = dpi,
width = 90,
height = 45,
units = "mm")Report MAE
- Average marginal effect -0.08 SDs (95 % CI [-0.11, -0.06])
Methods
Task and virtual environment
Trial types of the experiment
Encoding trial
Average trial lengths:
# Subset data to encoding only
encoding_only <- OLM_7T_trial_results[OLM_7T_trial_results$trialType == "encoding", ]
# Calculate trial duration
encoding_only$trial_duration <- encoding_only$end_time - encoding_only$start_time
# Calculate the average for each subject
trial_dur <- ddply(encoding_only, c("subject"), summarise,
total_length = sum(trial_duration),
run_length = total_length/2,
trial_duration = mean(trial_duration),
N = length(subject))
# Create report strings
str1 <- mean_SD_str2(trial_dur$trial_duration, report_type, digits1, rounding_type)
str2 <- mean_SD_str2(trial_dur$run_length, report_type, digits1, rounding_type)Average trial duration of an encoding trial: M = 20 (SD = 3.1).
Average encoding run duration: M = 360 (SD = 55).
Procedure
Days between sessions
# Load lab notebook with information on each subject
lab_notebook <- read_excel("data/ignore_fMRI_version1/VP00157_PPT_Info_Final.xlsx")
# Subset to only the participants used in this analysis
lab_notebook <- lab_notebook[lab_notebook$`ID Number` %in% subjects, ]
# Extract time between sessions
val <- lab_notebook$`Delay between sessions`
# Create report string
str1 <- mean_SD_str2(val, report_type, digits1, rounding_type)- Number of days between sessions M = 12 (SD = 13)
MRI acquisition
Number of volumes
# Load the .RData with the values
load("data/ignore_fMRI_version1/extracted_values/SpaNov_nVol_mean.RData")
# Calculate the total number of volumes
nVol_mean$nVol_total <- nVol_mean$nVol1 + nVol_mean$nVol2
# Convert this to minutes
nVol_mean$dataMinutes <- nVol_mean$nVol_total/60
# Create report strings
str1 <- mean_SD_str2(nVol_mean$nVol1, report_type, digits1, rounding_type)
str2 <- mean_SD_str2(nVol_mean$nVol2, report_type, digits1, rounding_type)
minTime <- signif(min(nVol_mean$dataMinutes), 3)- Number of volumes Run 1: M = 410 (SD = 81)
- Number of volumes Run 2: M = 370 (SD = 37)
- Minimum data in minutes: 10.5
Neuroimaging analysis overview
Standard GLMs
Gradient analysis
Hippocampal gradient analysis
Get the number of average voxels per position
voxelNum <- mean_SD_str2(HC_data_agg_pos$n, 1, digits1, rounding_type)The average number of voxels per position is M = 7.6 (SD = 3)
Supplementary material
Supplementary notes
Behavioural models of novelty separately
# Visits
load("fitted_brms_models/SpaNov_m_visits_poisson.Rdata")
mae1 <- avg_comparisons(m_visits_run1, variables = "s_onset_rel", type = "response", re_formula = NULL)
mae2 <- avg_comparisons(m_visits_run2, variables = "s_onset_rel", type = "response", re_formula = NULL)
# Time since
load("fitted_brms_models/SpaNov_m_timeSince.Rdata")
mae3 <- avg_comparisons(m_timeSince_run1, variables = "s_onset_rel", type = "response", re_formula = NULL)
mae4 <- avg_comparisons(m_timeSince_run2, variables = "s_onset_rel", type = "response", re_formula = NULL)
# Make strings
str1 <- create_str_from_avg_comparisons(mae1, "visits")
str2 <- create_str_from_avg_comparisons(mae2, "visits")
str3 <- create_str_from_avg_comparisons(mae3, "seconds")
str4 <- create_str_from_avg_comparisons(mae4, "seconds")Average marginal effects: - Run 1: 1.85 visits (95 % CI [1.63, 2.09]) - Run 2: 2.16 visits (95 % CI [1.95, 2.37]) - Run 1: 56.34 seconds (95 % CI [49.46, 63.84]) - Run 2: -145.42 seconds (95 % CI [-157.7, -133.73])
Supplementary tables
Significant clusters for linear contrast (Level 1 < Level 2 < Level 3 < …)
Click here for table.
| part | cluster | size | zValue_mean | zValue_median | zValue_max | zValue_min | zValue_peak_abs | cluster_mass | clsuter_mass_log | peak_x | peak_y | peak_z | peak_region | peak_region_name | num_regions |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| cortex_left | 1 | 22 | 3.063706 | 2.999546 | 3.6 | 2.5 | 3.638579 | 67 | 4.210668 | -3.791392 | -46.5213089 | 30.8302307 | 214 | d23ab | 2 |
| cortex_left | 2 | 479 | 3.210936 | 3.128089 | 4.6 | 2.5 | 4.618636 | 1500 | 7.338263 | -22.350677 | -43.2125931 | 72.1212463 | 232 | 2 | 7 |
| cortex_left | 3 | 357 | 3.129615 | 3.080705 | 4.4 | 2.5 | 4.356853 | 1100 | 7.018646 | -44.840401 | -25.9435463 | 17.5734921 | 284 | RI | 8 |
| cortex_left | 4 | 1425 | 3.630920 | 3.513689 | 6.4 | 2.5 | 6.428157 | 5200 | 8.551413 | -5.515844 | -16.9644909 | 58.8415260 | 220 | 24dd | 16 |
| cortex_left | 9 | 30 | 2.695159 | 2.697803 | 2.9 | 2.5 | 2.874365 | 81 | 4.392655 | -9.651464 | -42.0632935 | 60.6421967 | 217 | 5mv | 2 |
| cortex_left | 12 | 468 | 3.346280 | 3.275896 | 5.1 | 2.5 | 5.148257 | 1600 | 7.356318 | -51.256317 | -17.4799042 | 51.5853081 | 189 | 3b | 6 |
| cortex_left | 18 | 44 | 2.766905 | 2.739925 | 3.1 | 2.5 | 3.093434 | 120 | 4.801919 | -59.008034 | -37.4201393 | -0.5577022 | 309 | STSdp | 2 |
| cortex_left | 19 | 70 | 3.043194 | 2.963094 | 3.9 | 2.6 | 3.897602 | 210 | 5.361403 | -61.696888 | -43.7761002 | 11.0784178 | 319 | TPOJ1 | 2 |
| cortex_left | 20 | 28 | 2.984275 | 2.937808 | 3.5 | 2.5 | 3.458510 | 84 | 4.425561 | -61.709724 | -46.8277969 | -2.5455508 | 313 | TE1p | 2 |
| cortex_left | 24 | 46 | 3.100437 | 3.043691 | 3.7 | 2.5 | 3.705008 | 140 | 4.960185 | -44.560581 | -14.6814651 | 0.5779716 | 283 | 52 | 4 |
| cortex_left | 26 | 272 | 3.135538 | 2.986855 | 4.8 | 2.5 | 4.801907 | 850 | 6.748603 | -58.528763 | 9.8367786 | 16.2200336 | 258 | 6r | 8 |
| cortex_left | 39 | 109 | 3.257279 | 3.251825 | 4.2 | 2.6 | 4.240385 | 360 | 5.872240 | -50.125221 | -73.3484497 | 8.6767569 | 182 | MST | 4 |
| cortex_left | 40 | 180 | 3.216062 | 3.199615 | 4.4 | 2.5 | 4.369700 | 580 | 6.361115 | -51.717174 | -66.0832825 | 29.1244259 | 330 | PGi | 2 |
| cortex_left | 42 | 23 | 2.880656 | 2.844342 | 3.3 | 2.5 | 3.317583 | 66 | 4.193512 | -54.612988 | -0.3613958 | 42.1993637 | 236 | 6v | 4 |
| cortex_left | 43 | 77 | 2.885721 | 2.847735 | 3.4 | 2.5 | 3.436671 | 220 | 5.403580 | -62.776825 | -33.4391670 | 32.7153397 | 328 | PF | 3 |
| cortex_left | 47 | 121 | 2.972278 | 2.908775 | 3.8 | 2.4 | 3.818675 | 360 | 5.885119 | -7.869712 | 47.3418388 | -12.9690065 | 245 | 10r | 4 |
| cortex_left | 55 | 150 | 3.014078 | 2.915844 | 4.0 | 2.5 | 4.010159 | 450 | 6.113929 | -55.074410 | 4.9239883 | -23.9912205 | 308 | STSda | 5 |
| cortex_left | 65 | 279 | 3.013078 | 2.995985 | 3.8 | 2.4 | 3.759457 | 840 | 6.734174 | -13.645574 | 46.4050636 | 47.3006058 | 251 | 9p | 7 |
| cortex_left | 68 | 29 | 3.141467 | 3.073398 | 4.0 | 2.5 | 4.037845 | 91 | 4.511986 | -8.394770 | 51.6081810 | 7.1561594 | 249 | 9m | 3 |
| cortex_left | 71 | 42 | 2.948894 | 2.914134 | 3.6 | 2.5 | 3.602218 | 120 | 4.819100 | -34.522411 | 49.6495667 | 31.6823597 | 264 | 46 | 1 |
| cortex_right | 87 | 1643 | 3.457321 | 3.398987 | 5.7 | 2.4 | 5.702450 | 5700 | 8.644773 | 13.677155 | -16.6584091 | 75.3533707 | 55 | 6mp | 20 |
| cortex_right | 92 | 505 | 3.332931 | 3.249488 | 5.0 | 2.5 | 5.003071 | 1700 | 7.428411 | 33.398880 | -40.5589371 | 66.4627914 | 52 | 2 | 6 |
| cortex_right | 93 | 31 | 2.747490 | 2.681356 | 3.2 | 2.5 | 3.156602 | 85 | 4.444675 | 30.130425 | -26.7811832 | 67.4318771 | 53 | 3a | 2 |
| cortex_right | 95 | 20 | 2.676820 | 2.652724 | 2.9 | 2.5 | 2.918891 | 54 | 3.980362 | 22.743727 | -42.0650024 | 73.2298584 | 52 | 2 | 1 |
| cortex_right | 96 | 141 | 3.187271 | 3.139549 | 4.2 | 2.5 | 4.176615 | 450 | 6.107925 | 57.448917 | -13.8404722 | 48.2410164 | 51 | 1 | 3 |
| cortex_right | 99 | 457 | 3.089366 | 2.971921 | 4.8 | 2.4 | 4.844986 | 1400 | 7.252649 | 47.416145 | -25.7598839 | 18.3474407 | 104 | RI | 9 |
| cortex_right | 102 | 80 | 3.276918 | 3.268695 | 4.5 | 2.6 | 4.506863 | 260 | 5.568930 | 62.873524 | -40.0575027 | 15.5642920 | 28 | STV | 3 |
| cortex_right | 103 | 37 | 2.822053 | 2.749717 | 3.4 | 2.4 | 3.436898 | 100 | 4.648383 | 57.842972 | -34.0697250 | -1.2721851 | 129 | STSdp | 2 |
| cortex_right | 110 | 21 | 2.753690 | 2.766122 | 3.2 | 2.5 | 3.189069 | 58 | 4.057464 | 46.799515 | -2.4484978 | -9.9124584 | 167 | PoI1 | 2 |
| cortex_right | 112 | 55 | 3.004002 | 2.978083 | 3.7 | 2.6 | 3.667743 | 170 | 5.107278 | 45.115521 | 8.0075855 | 8.4428720 | 115 | FOP2 | 2 |
| cortex_right | 113 | 25 | 2.951910 | 2.927984 | 3.4 | 2.5 | 3.383995 | 74 | 4.301328 | 41.860241 | 29.4441929 | 3.9625542 | 108 | FOP4 | 2 |
| cortex_right | 116 | 37 | 3.206650 | 3.216331 | 3.9 | 2.4 | 3.909213 | 120 | 4.776145 | 16.628929 | -83.4581909 | 40.5141487 | 3 | V6 | 2 |
| cortex_right | 117 | 220 | 3.506168 | 3.399297 | 5.7 | 2.4 | 5.696318 | 770 | 6.648151 | 50.564453 | -72.4114990 | 6.8142180 | 23 | MT | 7 |
| cortex_right | 125 | 216 | 3.174723 | 3.072107 | 4.3 | 2.5 | 4.313631 | 690 | 6.530499 | 61.787083 | 8.8548021 | 19.2581615 | 56 | 6v | 6 |
| cortex_right | 130 | 24 | 2.946266 | 2.877828 | 3.7 | 2.5 | 3.672417 | 71 | 4.258592 | 43.409565 | 17.7264061 | -35.5937271 | 131 | TGd | 1 |
| cortex_right | 146 | 38 | 2.828438 | 2.834183 | 3.3 | 2.5 | 3.254191 | 110 | 4.677311 | 8.638423 | 55.7692184 | 25.6142464 | 69 | 9m | 2 |
| cortex_right | 147 | 25 | 2.989004 | 2.958746 | 3.6 | 2.5 | 3.600116 | 75 | 4.313816 | 35.940197 | 50.5307732 | 30.8277359 | 84 | 46 | 1 |
| cortex_right | 155 | 26 | 2.818180 | 2.809673 | 3.1 | 2.6 | 3.136499 | 73 | 4.294188 | 58.314800 | 0.0648438 | -16.3995800 | 128 | STSda | 1 |
| subcort | 168 | 161 | 3.074034 | 2.976389 | 4.5 | 2.4 | 4.542136 | 490 | 6.204395 | 26.000000 | -88.0000000 | -38.0000000 | 371 | CEREBELLUM_RIGHT | 1 |
| subcort | 193 | 107 | 2.998622 | 2.911509 | 4.1 | 2.5 | 4.123521 | 320 | 5.770982 | 20.000000 | -14.0000000 | -20.0000000 | 376 | HIPPOCAMPUS_RIGHT | 2 |
| subcort | 194 | 135 | 3.144796 | 3.095737 | 4.5 | 2.5 | 4.509548 | 420 | 6.051024 | -18.000000 | -12.0000000 | -20.0000000 | 367 | HIPPOCAMPUS_LEFT | 2 |
| subcort | 215 | 43 | 2.861746 | 2.787239 | 3.8 | 2.4 | 3.830762 | 120 | 4.812632 | -16.000000 | 8.0000000 | -8.0000000 | 364 | PUTAMEN_LEFT | 1 |
| part | cluster | regions | labels |
|---|---|---|---|
| cortex_left | 1 | 215,214 | 31pv,d23ab |
| cortex_left | 2 | 232,231,227,222,229,219,225 | 2,1,7PC,7AL,VIP,5L,7Am |
| cortex_left | 3 | 282,328,285,284,204,205,354,281 | OP2-3,PF,PFcm,RI,A1,PSL,LBelt,OP1 |
| cortex_left | 4 | 218,221,223,240,237,239,189,233,188,235,224,216,217,220,234,276 | 23c,24dv,SCEF,p32pr,p24pr,a24pr,3b,3a,4,6mp,6ma,5m,5mv,24dd,6d,6a |
| cortex_left | 9 | 219,217 | 5L,5mv |
| cortex_left | 12 | 189,233,188,190,234,231 | 3b,3a,4,FEF,6d,1 |
| cortex_left | 18 | 310,309 | STSvp,STSdp |
| cortex_left | 19 | 319,208 | TPOJ1,STV |
| cortex_left | 20 | 310,313 | STSvp,TE1p |
| cortex_left | 24 | 282,283,347,348 | OP2-3,52,PoI1,Ig |
| cortex_left | 26 | 288,295,294,293,258,279,236,188 | FOP4,FOP2,FOP3,FOP1,6r,43,6v,4 |
| cortex_left | 39 | 203,182,320,321 | MT,MST,TPOJ2,TPOJ3 |
| cortex_left | 40 | 330,331 | PGi,PGs |
| cortex_left | 42 | 188,236,191,192 | 4,6v,PEF,55b |
| cortex_left | 43 | 328,327,285 | PF,PFop,PFcm |
| cortex_left | 47 | 268,245,345,252 | 10v,10r,s32,10d |
| cortex_left | 55 | 311,303,308,356,312 | TGd,STGa,STSda,STSva,TE1a |
| cortex_left | 65 | 206,250,251,252,249,267,248 | SFL,8BL,9p,10d,9m,9a,8Ad |
| cortex_left | 68 | 241,244,249 | a24,p32,9m |
| cortex_left | 71 | 264 | 46 |
| cortex_right | 87 | 12,38,41,43,60,57,8,55,44,37,39,36,40,52,51,10,9,54,53,96 | 55b,23c,24dv,SCEF,p32pr,p24pr,4,6mp,6ma,5mv,5ROI,5m,24dd,2,1,FEF,3b,6d,3a,6a |
| cortex_right | 92 | 39,52,47,42,49,48 | 5ROI,2,7PC,7AROI,VIP,LIPv |
| cortex_right | 93 | 9,53 | 3b,3a |
| cortex_right | 95 | 52 | 2 |
| cortex_right | 96 | 51,52,116 | 1,2,PFt |
| cortex_right | 99 | 148,105,104,25,124,174,101,147,100 | PF,PFcm,RI,PSROI,PBelt,LBelt,OP1,PFop,OP4 |
| cortex_right | 102 | 139,28,25 | TPOJ1,STV,PSROI |
| cortex_right | 103 | 129,125 | STSdp,A5 |
| cortex_right | 110 | 167,178 | PoI1,PI |
| cortex_right | 112 | 115,114 | FOP2,FOP3 |
| cortex_right | 113 | 108,169 | FOP4,FOP5 |
| cortex_right | 116 | 152,3 | V6A,V6 |
| cortex_right | 117 | 156,2,23,137,157,140,141 | V4t,MST,MT,PHT,FST,TPOJ2,TPOJ3 |
| cortex_right | 125 | 113,78,99,56,11,12 | FOP1,6r,43,6v,PEF,55b |
| cortex_right | 130 | 131 | TGd |
| cortex_right | 146 | 87,69 | 9a,9m |
| cortex_right | 147 | 84 | 46 |
| cortex_right | 155 | 128 | STSda |
| subcort | 168 | 371 | CEREBELLUM_RIGHT |
| subcort | 193 | 376,377 | HIPPOCAMPUS_RIGHT,AMYGDALA_RIGHT |
| subcort | 194 | 367,368 | HIPPOCAMPUS_LEFT,AMYGDALA_LEFT |
| subcort | 215 | 364 | PUTAMEN_LEFT |
Significant clusters for linear contrast (Level 1 > Level 2 > Level 3 > …)
Click here for table.
| part | cluster | size | zValue_mean | zValue_median | zValue_max | zValue_min | zValue_peak_abs | cluster_mass | clsuter_mass_log | peak_x | peak_y | peak_z | peak_region | peak_region_name | num_regions |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| cortex_left | 1 | 90 | -3.519341 | -3.443073 | -2.4 | -4.8 | 4.790668 | 320 | 5.758084 | -2.259454 | -31.05509 | 32.052841 | 194 | RSC | 3 |
| cortex_left | 2 | 2060 | -3.930611 | -3.876594 | -2.5 | -6.7 | 6.735095 | 8100 | 8.999256 | -32.472549 | -88.98758 | 29.779585 | 326 | IP0 | 28 |
| cortex_left | 3 | 221 | -3.194249 | -3.094511 | -2.5 | -4.3 | 4.344084 | 710 | 6.559515 | -28.574898 | -52.22017 | -18.019115 | 307 | PHA3 | 10 |
| cortex_left | 5 | 27 | -2.856392 | -2.848945 | -2.5 | -3.3 | 3.321889 | 77 | 4.345396 | -8.653727 | -54.54932 | 5.236495 | 211 | POS1 | 2 |
| cortex_left | 7 | 84 | -3.043542 | -3.026831 | -2.5 | -3.7 | 3.736013 | 260 | 5.543839 | -5.988591 | 25.43487 | 52.024731 | 243 | 8BM | 4 |
| cortex_left | 11 | 216 | -3.485387 | -3.370426 | -2.5 | -5.1 | 5.099926 | 750 | 6.623857 | -27.186953 | 19.61117 | 58.827522 | 277 | i6-8 | 8 |
| cortex_left | 14 | 107 | -3.212653 | -3.204306 | -2.5 | -4.1 | 4.055283 | 340 | 5.839926 | -47.295277 | 16.06033 | 34.457317 | 260 | IFJp | 4 |
| cortex_left | 15 | 131 | -3.085281 | -3.057481 | -2.6 | -3.9 | 3.863326 | 400 | 6.001840 | -44.775909 | 41.61255 | 25.340578 | 263 | p9-46v | 4 |
| cortex_left | 17 | 55 | -3.049246 | -2.937280 | -2.5 | -3.8 | 3.841370 | 170 | 5.122228 | -36.455795 | 60.54835 | 7.798744 | 265 | a9-46v | 3 |
| cortex_left | 29 | 483 | -4.251315 | -4.164235 | -2.5 | -6.8 | 6.760676 | 2100 | 7.627245 | -13.782526 | -92.53966 | -12.004962 | 184 | V2 | 4 |
| cortex_left | 32 | 22 | -2.944452 | -2.925727 | -2.5 | -3.4 | 3.412598 | 65 | 4.170965 | -6.316596 | -90.02340 | 8.658187 | 181 | V1 | 1 |
| cortex_right | 38 | 2303 | -3.929091 | -3.748647 | -2.4 | -6.9 | 6.878045 | 9000 | 9.110376 | 36.958008 | -87.18770 | 24.761272 | 146 | IP0 | 35 |
| cortex_right | 39 | 244 | -3.485530 | -3.477159 | -2.5 | -4.9 | 4.867373 | 850 | 6.745788 | 24.893423 | -51.38643 | -16.057383 | 160 | VMV2 | 7 |
| cortex_right | 40 | 118 | -3.852320 | -3.901859 | -2.5 | -4.9 | 4.903612 | 450 | 6.119360 | 3.001509 | -31.90833 | 32.272301 | 14 | RSC | 1 |
| cortex_right | 44 | 63 | -3.120948 | -3.171131 | -2.5 | -3.8 | 3.839579 | 200 | 5.281271 | 6.831443 | 29.15144 | 44.001465 | 179 | a32pr | 3 |
| cortex_right | 46 | 402 | -3.335551 | -3.302893 | -2.5 | -4.7 | 4.685742 | 1300 | 7.201090 | 31.845680 | 21.57755 | 58.013378 | 97 | i6-8 | 10 |
| cortex_right | 50 | 186 | -3.224232 | -3.150901 | -2.5 | -4.5 | 4.506076 | 600 | 6.396442 | 48.756924 | 39.91827 | 20.827803 | 81 | IFSp | 4 |
| cortex_right | 70 | 472 | -4.438338 | -4.449830 | -2.4 | -7.2 | 7.179883 | 2100 | 7.647259 | 13.675855 | -91.91169 | -6.089150 | 1 | V1 | 5 |
| cortex_right | 74 | 28 | -2.969049 | -2.899777 | -2.6 | -3.6 | 3.633223 | 83 | 4.420446 | 10.256091 | -91.39706 | 9.405240 | 1 | V1 | 2 |
| cortex_right | 77 | 58 | -3.031525 | -2.926241 | -2.5 | -4.5 | 4.476549 | 180 | 5.169509 | 23.169453 | 66.53020 | 2.306681 | 170 | p10p | 3 |
| cortex_right | 80 | 23 | -2.992635 | -3.061711 | -2.5 | -3.5 | 3.546734 | 69 | 4.231649 | 26.927168 | 42.81836 | 44.168594 | 68 | 8Ad | 2 |
| subcort | 82 | 172 | -3.229611 | -3.108652 | -2.5 | -4.8 | 4.768074 | 560 | 6.319856 | -38.000000 | -70.00000 | -52.000000 | 361 | CEREBELLUM_LEFT | 1 |
| subcort | 84 | 283 | -3.336223 | -3.222986 | -2.4 | -5.3 | 5.324766 | 940 | 6.850286 | 10.000000 | -46.00000 | -48.000000 | 371 | CEREBELLUM_RIGHT | 2 |
| subcort | 86 | 185 | -3.337890 | -3.209643 | -2.4 | -4.7 | 4.740293 | 620 | 6.425695 | 34.000000 | -70.00000 | -50.000000 | 371 | CEREBELLUM_RIGHT | 1 |
| subcort | 87 | 198 | -3.299109 | -3.250997 | -2.5 | -4.8 | 4.785460 | 650 | 6.481920 | -24.000000 | -32.00000 | -42.000000 | 361 | CEREBELLUM_LEFT | 1 |
| subcort | 92 | 372 | -3.564598 | -3.446843 | -2.5 | -5.7 | 5.698634 | 1300 | 7.189945 | -4.000000 | -80.00000 | -40.000000 | 361 | CEREBELLUM_LEFT | 1 |
| subcort | 99 | 271 | -3.506712 | -3.379627 | -2.5 | -6.0 | 6.013061 | 950 | 6.856798 | 8.000000 | -78.00000 | -24.000000 | 371 | CEREBELLUM_RIGHT | 1 |
| subcort | 106 | 20 | -3.091158 | -3.096199 | -2.6 | -4.0 | 4.030955 | 62 | 4.124278 | -10.000000 | -58.00000 | -38.000000 | 361 | CEREBELLUM_LEFT | 1 |
| subcort | 112 | 201 | -3.193322 | -3.150749 | -2.5 | -4.6 | 4.587187 | 640 | 6.464367 | 28.000000 | -64.00000 | -30.000000 | 371 | CEREBELLUM_RIGHT | 1 |
| subcort | 113 | 250 | -3.332495 | -3.219638 | -2.5 | -5.0 | 4.964678 | 830 | 6.725182 | -26.000000 | -60.00000 | -32.000000 | 361 | CEREBELLUM_LEFT | 1 |
| subcort | 141 | 47 | -2.979207 | -2.928996 | -2.5 | -4.0 | 4.010671 | 140 | 4.941805 | 18.000000 | -32.00000 | 6.000000 | 372 | THALAMUS_RIGHT | 1 |
| part | cluster | regions | labels |
|---|---|---|---|
| cortex_left | 1 | 194,214,212 | RSC,d23ab,23d |
| cortex_left | 2 | 211,301,207,218,342,297,324,329,230,197,195,322,226,225,209,210,199,326,323,331,325,275,228,186,338,321,330,184 | POS1,ProS,PCV,23c,31a,AIP,IP2,PFm,MIP,IPS1,POS2,DVT,7PL,7Am,7Pm,7m,V3B,IP0,PGp,PGs,IP1,LIPd,LIPv,V4,V3CD,TPOJ3,PGi,V2 |
| cortex_left | 3 | 306,333,301,299,184,340,307,335,343,334 | PHA1,VMV1,ProS,PreS,V2,VMV2,PHA3,PHA2,VVC,VMV3 |
| cortex_left | 5 | 194,211 | RSC,POS1 |
| cortex_left | 7 | 206,243,223,359 | SFL,8BM,SCEF,a32pr |
| cortex_left | 11 | 190,276,224,278,248,277,247,192 | FEF,6a,6ma,s6-8,8Ad,i6-8,8Av,55b |
| cortex_left | 14 | 253,259,260,191 | 8C,IFJa,IFJp,PEF |
| cortex_left | 15 | 253,263,262,261 | 8C,p9-46v,IFSa,IFSp |
| cortex_left | 17 | 265,266,350 | a9-46v,9-46d,p10p |
| cortex_left | 29 | 181,184,185,186 | V1,V2,V3,V4 |
| cortex_left | 32 | 181 | V1 |
| cortex_right | 38 | 149,31,121,14,162,27,35,38,34,144,117,50,17,30,15,142,46,45,152,29,161,6,19,146,143,151,13,145,95,48,20,158,150,1,3 | PFm,POS1,ProS,RSC,31a,PCV,31pv,23c,d23ab,IP2,AIP,MIP,IPS1,7m,POS2,DVT,7PROI,7Am,V6A,7Pm,31pd,V4,V3B,IP0,PGp,PGs,V3A,IP1,LIPd,LIPv,LO1,V3CD,PGi,V1,V6 |
| cortex_right | 39 | 160,126,153,127,155,163,154 | VMV2,PHA1,VMV1,PHA3,PHA2,VVC,VMV3 |
| cortex_right | 40 | 14 | RSC |
| cortex_right | 44 | 63,43,179 | 8BM,SCEF,a32pr |
| cortex_right | 46 | 10,96,73,80,79,11,98,68,97,67 | FEF,6a,8C,IFJp,IFJa,PEF,s6-8,8Ad,i6-8,8Av |
| cortex_right | 50 | 73,83,81,82 | 8C,p9-46v,IFSp,IFSa |
| cortex_right | 70 | 1,4,5,6,7 | V1,V2,V3,V4,V8 |
| cortex_right | 74 | 1,4 | V1,V2 |
| cortex_right | 77 | 170,90,72 | p10p,10pp,10d |
| cortex_right | 80 | 68,84 | 8Ad,46 |
| subcort | 82 | 361 | CEREBELLUM_LEFT |
| subcort | 84 | 371,366 | CEREBELLUM_RIGHT,BRAIN_STEM |
| subcort | 86 | 371 | CEREBELLUM_RIGHT |
| subcort | 87 | 361 | CEREBELLUM_LEFT |
| subcort | 92 | 361 | CEREBELLUM_LEFT |
| subcort | 99 | 371 | CEREBELLUM_RIGHT |
| subcort | 106 | 361 | CEREBELLUM_LEFT |
| subcort | 112 | 371 | CEREBELLUM_RIGHT |
| subcort | 113 | 361 | CEREBELLUM_LEFT |
| subcort | 141 | 372 | THALAMUS_RIGHT |
Supplementary figures
Learning curves
# Creating these figures takes too long so eval = FALSE
load("fitted_brms_models/SpaNov_m_navTime.Rdata")
# Subset to retrieval trials
OLM_7T_encoding <- OLM_7T_trial_results[OLM_7T_trial_results$trialType == "encoding", ]
# Calculate distance between start and object
x1 <- OLM_7T_encoding$start_x
z1 <- OLM_7T_encoding$start_z
x2 <- OLM_7T_encoding$object_x
z2 <- OLM_7T_encoding$object_z
OLM_7T_encoding$start2obj <- euclideanDistance3D(x1, 1, z1, x2, 1, z2)
# Scale values
OLM_7T_encoding$s_timesObjectPresented <- scale_m0_sd0.5(OLM_7T_encoding$timesObjectPresented)
OLM_7T_encoding$s_start2obj <- scale_m0_sd0.5(OLM_7T_encoding$start2obj)
# Rename model & set variable of interest
m <- m_navTime3
voi <- "s_timesObjectPresented"
# Get which raw values the scaled values correspond to
raw_and_scaled_values <- ddply(OLM_7T_encoding, c("timesObjectPresented"), summarise,
raw = timesObjectPresented[1],
scaled = mean(s_timesObjectPresented))
# Calculate marginal average effect is the contrast 0 vs. 1
mae1_navTime_respone <- avg_comparisons(m, variables = voi, type = "response", re_formula = NULL)
#mae1_navTime_link <- avg_comparisons(m, variables = voi, type = "link", re_formula = NULL)
# Get posterior draws
mae1_navTime_respone_draws <- get_draws(mae1_navTime_respone)
p1_mae_draws <- ggplot(mae1_navTime_respone_draws, aes(x = draw)) +
geom_vline(xintercept = 0, linetype = "dashed") +
stat_halfeye(fill = base_col[1], point_size = 1.5) +
labs(x = "Navigation time (seconds)", y = "Density",
title = "Marginal average effect\n of number of times presented") +
theme_journal() +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
# Create estimates for line
grid <- datagrid(s_timesObjectPresented = seq_min_2_max(raw_and_scaled_values$scaled),
subject = unique, s_start2obj = 0, model = m)
pred_line_navTime <- predictions(m, newdata = grid, by = voi, re_formula = NULL, type = "response")
pred_line_navTime <- as.data.frame(pred_line_navTime)
# Create estimates actually existing values
grid <- datagrid(s_timesObjectPresented = unique, subject = unique, s_start2obj = 0, model = m)
pred_actual_navTime <- predictions(m, newdata = grid, by = voi, re_formula = NULL, type = "response")
pred_actual_navTime <- as.data.frame(pred_actual_navTime)
# Calculate subject-level lines
grid <- datagrid(s_timesObjectPresented = seq_min_2_max(raw_and_scaled_values$scaled),
subject = unique, s_start2obj = 0, model = m)
pred_line_subject_navTime <- predictions(m, newdata = grid, type = "response")
pred_line_subject_navTime <- as.data.frame(pred_line_subject_navTime)
# Make DFs unique
pred_line_navTime1 <- pred_line_navTime
pred_actual_navTime1 <- pred_actual_navTime
pred_line_subject_navTime1 <- pred_line_subject_navTime
# Assemble plot
p1_line <-ggplot() +
geom_line(data = pred_line_subject_navTime1, mapping = aes(x = s_timesObjectPresented, y = estimate, group = subject),
colour = "darkgrey") +
geom_line(data = pred_line_navTime1, mapping = aes(x = s_timesObjectPresented, y = estimate), linewidth = 0.3) +
geom_point(data = pred_actual_navTime1, mapping = aes(x = s_timesObjectPresented, y = estimate)) +
geom_errorbar(data = pred_actual_navTime1, mapping = aes(x = s_timesObjectPresented, ymin = conf.low, ymax = conf.high),
width = 0.02) +
scale_x_continuous(breaks = raw_and_scaled_values$scaled, labels = raw_and_scaled_values$raw) +
theme_journal() +
labs(x = "Number of times the object was presented", y = "Navigation time\n(seconds)")
###########################################################################_
# Assemble figure
learning_curves <- plot_grid(p1_mae_draws, p1_line, align = "hv", ncol = 2,
labels = "auto", label_size = 10)
# Save
ggsave(learning_curves,
filename = paste0(figurePath, "Fig_S01_learning_curves.png"),
dpi = dpi, width = 130, height = 40, units = "mm")
Report MAEs
# Load the model
load("fitted_brms_models/SpaNov_m_navTime.Rdata")
# Rename model & set variable of interest
m <- m_navTime3
voi <- "s_timesObjectPresented"
mae1_navTime_respone <- avg_comparisons(m, variables = voi, type = "response", re_formula = NULL)
str1 <- create_str_from_avg_comparisons(mae1_navTime_respone, "seconds")- Average marginal effect:
str 1
Locomotion by novelty
# Load the image from the gradient event file creation
load("event_tables/images/SpaNov_event_file_gradients.RData")
# Create a data frame
ev_info <- rbindlist(condition_info)
# Subset to the analysis with 6 levels and encoding and for the novelty score
ev_info <- ev_info[ev_info$curr_numLvl == 6 & ev_info$runType == "e" & ev_info$type == "noveltyScore", ]
# Calculate average across runs
ev_info_agg <- ddply(ev_info, c("subj", "level"), summarise,
mean_per_tra_arc = mean(mean_per_tra_arc),
mean_per_rot_arc = mean(mean_per_rot_arc),
mean_per_sta_arc = mean(mean_per_sta_arc))
# Calculate average for each run
ev_info_agg1 <- ddply(ev_info, c("subj", "run"), summarise,
mean_per_tra_arc = mean(mean_per_tra_arc),
mean_per_rot_arc = mean(mean_per_rot_arc),
mean_per_sta_arc = mean(mean_per_sta_arc))
# Remove middle levels
ev_info_agg <- ev_info_agg[ev_info_agg$level == "lvl1" | ev_info_agg$level == "lvl6", ]
# Create the plots
## Translation
p1 <- ggplot(ev_info_agg, aes(x = level, y = mean_per_tra_arc, fill = level)) +
geom_line(aes(group = subj), colour = boxplot_line_col) +
geom_point(colour = boxplot_point_col) +
geom_boxplot(colour = boxplot_border_col, outlier.shape = NA, width = 0.4) +
stat_summary(geom = "point", fun = "mean", col = 'white', shape = 24,
position = position_dodge(width = 0.75), size = 0.8) +
scale_fill_manual(values = novFam_gradient[c(1, 6)]) +
theme_journal() +
scale_y_continuous(breaks = c(0.25, 0.82, 1.40)) +
coord_cartesian(xlim = c(0.5, 2.5), ylim = c(0.25, 1.40), expand = FALSE) +
theme(legend.position = "none") +
labs(title = "Translation", x = "Novelty level", y = "arcsine(percentage)")
## Rotation
p2 <- ggplot(ev_info_agg, aes(x = level, y = mean_per_rot_arc, fill = level)) +
geom_line(aes(group = subj), colour = boxplot_line_col) +
geom_point(colour = boxplot_point_col) +
geom_boxplot(colour = boxplot_border_col, outlier.shape = NA, width = 0.4) +
stat_summary(geom = "point", fun = "mean", col = 'white', shape = 24,
position = position_dodge(width = 0.75), size = 0.8) +
scale_fill_manual(values = novFam_gradient[c(1, 6)]) +
theme_journal() +
scale_y_continuous(breaks = c(-1.5, -1, -0.58)) +
coord_cartesian(xlim = c(0.5, 2.5), ylim = c(-1.5, -0.58), expand = FALSE) +
theme(legend.position = "none") +
labs(title = "Rotation", x = "Novelty level", y = "arcsine(percentage)")
## Stationary
p3 <- ggplot(ev_info_agg, aes(x = level, y = mean_per_sta_arc, fill = level)) +
geom_line(aes(group = subj), colour = boxplot_line_col) +
geom_point(colour = boxplot_point_col) +
geom_boxplot(colour = boxplot_border_col, outlier.shape = NA, width = 0.4) +
stat_summary(geom = "point", fun = "mean", col = 'white', shape = 24,
position = position_dodge(width = 0.75), size = 0.8) +
scale_fill_manual(values = novFam_gradient[c(1, 6)]) +
theme_journal() +
scale_y_continuous(breaks = c(-1.5, -1.10, -0.73)) +
coord_cartesian(xlim = c(0.5, 2.5), ylim = c(-1.5, -0.73), expand = FALSE) +
theme(legend.position = "none") +
labs(title = "Stationary", x = "Novelty level", y = "arcsine(percentage)")
# Combine and save
p_comb <- plot_grid(p1, p2, p3, ncol = 3, align = "hv")
ggsave(p_comb,
filename = paste0(figurePath, "Fig_S02_locomotion.png"),
dpi = dpi,
width = 90,
height = 40,
units = "mm")
# Calculate the test statistics
## Translation
val1 <- ev_info_agg[ev_info_agg$level == 'lvl1', "mean_per_tra_arc"]
val2 <- ev_info_agg[ev_info_agg$level == 'lvl6', "mean_per_tra_arc"]
diff_val <- val1 - val2
BF1 <- signif(reportBF(ttestBF(diff_val)), digits1)
d1 <- signif(mean(diff_val)/sd(diff_val), digits1)
## Rotation
val1 <- ev_info_agg[ev_info_agg$level == 'lvl1', "mean_per_rot_arc"]
val2 <- ev_info_agg[ev_info_agg$level == 'lvl6', "mean_per_rot_arc"]
diff_val <- val1 - val2
BF2 <- signif(reportBF(ttestBF(diff_val)), digits1)
d2 <- signif(mean(diff_val)/sd(diff_val), digits1)
## Stationary
val1 <- ev_info_agg[ev_info_agg$level == 'lvl1', "mean_per_sta_arc"]
val2 <- ev_info_agg[ev_info_agg$level == 'lvl6', "mean_per_sta_arc"]
diff_val <- val1 - val2
BF3 <- signif(reportBF(ttestBF(diff_val)), digits1)
d3 <- signif(mean(diff_val)/sd(diff_val), digits1)- Translation Level 1 vs. Level 6: BF10 = 3.5^{15}, d = -1.7
- Rotation Level 1 vs. Level 6: BF10 = 7.3^{16}, d = 1.9
- Stationary Level 1 vs. Level 6: BF10 = 9.4^{8}, d = 1.1
Spatial distribution of novelty scores
# Circle definitions
## Area of the arena
arenaDiameter <- 180
radius1 <- arenaDiameter/2
areaCircle1 <- pi*(radius1^2)
## Area of circle with half the radius
radius2 <- radius1/2
areaCircle2 <- pi*radius2^2
# Draw inner and outer circle
## Outer circle
circle1 <- circleFun(c(0, 0), arenaDiameter, npoints = 1000)
circle1$z <- circle1$y # Rename to avoid error that columns do not match
circle1$y <- NULL
## Inner circle
circle2 <- circleFun(c(0, 0), radius2*2, npoints = 1000)
circle2$z <- circle2$y # Rename to avoid error that columns do not match
circle2$y <- NULL
# Convert list to data frame
event_info_df <- rbindlist(event_info, idcol = "id")
event_info_df <- as.data.frame(event_info_df)
# Subset to gradient level 1
event_info_df_lvl1 <- event_info_df[event_info_df$noveltyScore_gradientLevel == "lvl1" &
!is.na(event_info_df$noveltyScore_gradientLevel), ]
event_info_df_lvl1 <- event_info_df_lvl1[event_info_df_lvl1$run %in% c(1, 3), ]
spa_nov_dist <- ggplot() +
# Circles
geom_polygon(data = circle1, mapping = aes(x = x, y = z), fill = "lightgrey") +
geom_polygon(data = circle2, mapping = aes(x = x, y = z), fill = "darkgrey") +
# Lines and text for outer radius
geom_segment(aes(x = 0, y = -97, xend = radius1, yend = -97)) +
geom_segment(aes(x = radius1, y = 0, xend = radius1, yend = -97), linetype = "dashed") +
geom_segment(aes(x = 0, y = 0, xend = 0, yend = -97), linetype = "dashed") +
annotate("text", x = 45, y = -102, label = c("90 vm"), size = 2) +
# Lines for inner radius
geom_segment(aes(x = 0, y = -95, xend = radius2, yend = -95)) +
geom_segment(aes(x = radius2, y = 0, xend = radius2, yend = -95), linetype = "dashed") +
annotate("text", x = radius2/2, y = -90, label = c("45 vm"), size = 2) +
# Points
geom_point(data = event_info_df_lvl1, mapping = aes(x = pos_x, y = pos_z)) +
# Arrows and text
geom_curve(aes(x = -95, y = 70, xend = -50, yend = 50), ncp = 10, arrow = arrow(length = unit(1, "mm"))) +
geom_curve(aes(x = 95, y = 70, xend = 30, yend = 0), ncp = 10, arrow = arrow(length = unit(1, "mm")), curvature = -0.5) +
annotate("text", x = c(-95, 95), y = c(73, 73), label = c("Periphery", "Centre"), size = 2) +
# Other settings
coord_equal(xlim = c(-105, 105)) +
theme_void() +
theme(plot.background = element_rect(fill = "white", colour = "white"))
ggsave(spa_nov_dist,
filename = paste0(figurePath, "Fig_S03_spa_nov_dist.png"),
dpi = dpi, width = 80, height = 80, units = "mm")
# Calculate centrality for each level
## Subset to encoding runs & none NA values
event_info_df_encoding <- event_info_df[event_info_df$run %in% c(1, 3), ]
event_info_df_encoding <- na.omit(event_info_df_encoding)
## Calculate distance to centre
event_info_df_encoding$dist2centre <- euclideanDistance3D(0, 1, 0, event_info_df_encoding$pos_x, 1, event_info_df_encoding$pos_z)
## Categorise centrality
event_info_df_encoding$centrality1 <- ifelse(event_info_df_encoding$dist2centre >= radius2, 'perihery', 'central')
## Calculate average for each level
avg_level_centrality <- ddply(event_info_df_encoding, c("noveltyScore_gradientLevel"), summarise,
centrality1_score = mean(centrality1 == "central"))
Minimal influence of centrality on spatial novelty effects
# Loading files
dist2centre_file <- "/media/alex/work/Seafile/imaging_results/SpaNov/OLMe_7T_SpaNov_gradient_dist2centre-corrected_6lvl_smo4_MSMAll/cope7.feat/stats/vwc/results_lvl2cope1_dat_ztstat_c1.dscalar.nii"
SpaNov1_file <- "/media/alex/work/Seafile/imaging_results/SpaNov/OLMe_7T_SpaNov_gradient_6lvl_cue-delay_smo4_MSMAll/cope7.feat/stats/vwc/results_lvl2cope1_dat_ztstat_c1.dscalar.nii"
dist2centre_xii <- read_cifti(dist2centre_file)
SpaNov1_xii <- read_cifti(SpaNov1_file)
#
df <- data.frame(dist2centre = as.matrix(dist2centre_xii),
SpaNov1 = as.matrix(SpaNov1_xii))
# Create correlation string
r_val <- round(cor.test(df$dist2centre, df$SpaNov1)$estimate, 2)
p_val <- pValue(cor.test(df$dist2centre, df$SpaNov1)$p.value)
report_str <- paste0("r = ", r_val, ", p ", p_val)
# Create scatter plot
p_map_corr <- ggplot(df, aes(x = dist2centre, y = SpaNov1)) +
geom_density_2d_filled() +
#geom_smooth(method = "lm", formula = "y ~ x", colour = "white", se = FALSE) +
theme_journal() +
theme(legend.position = "none") +
theme(plot.margin = unit(c(1, 4, 1, 1), "mm")) +
coord_cartesian(expand = FALSE) +
annotate("text", x = I(0.3), y = I(0.9), label = report_str, size = 2, colour = "white") +
labs(x = "Z-map corrected for centrality", y = "Original z-map")
##############################################################################_
# Calculate centrality
## Area of the arena
radius1 <- 180/2
areaCircle1 <- pi*(radius1^2)
## Area of circle with half the radius
radius2 <- radius1/2
areaCircle2 <- pi*radius2^2
## Subset position data only encoding (include cue & delay because they are the same for all)
centrality_data <- locomotion_data
## Calculate distance to centre
centrality_data$dist2centre <- euclideanDistance3D(0, 1, 0, centrality_data$pos_x, 1, centrality_data$pos_z)
## Categorise centrality
centrality_data$centrality1 <- ifelse(centrality_data$dist2centre >= radius2, 'perihery', 'central')
## Calculate average
avg_centrality <- ddply(centrality_data, c("ppid","subject"), summarise,
score1 = mean(centrality1 == "central"))
## Write .csv for PALM analysis
write.csv(avg_centrality[,-2], file = "data/ignore_fMRI_version1/SpaNov_centrality_designMatrix_data.csv",
row.names = FALSE, quote = FALSE)
##############################################################################_
p <- plot_grid(p_map_corr, NULL, NULL,
nrow = 1, labels = "auto", label_size = 10, align = "hv")
# Combine plots
ggsave(p,
filename = paste0(figurePath, "Fig_S04_control_for_centrality.png"),
dpi = dpi, width = 180, height = 57, units = "mm")Individual differences in centrality do not predict spatial novelty
## Calculate average centrality
str1 <- mean_SD_str2(avg_centrality$score1, report_type, digits1, rounding_type)
## Create plot
avg_centrality_plot <- ggplot(avg_centrality, aes(x = score1)) +
geom_histogram(colour = base_col[1], fill = base_col[1]) +
geom_vline(xintercept = mean(avg_centrality$score1), colour = "black", linetype = 2, size = 0.5) +
scale_x_continuous(breaks = c(0.45, 0.55, 0.65)) +
scale_y_continuous(breaks = c(0, 5, 10)) +
coord_cartesian(xlim = c(0.44, 0.65), ylim = c(0, 10),
expand = FALSE) +
labs(x = "Centrality score", y = "Count") +
theme_journal() +
theme(plot.margin = unit(c(1, 2.5, 1, 2.5), "mm"))
##############################################################################_
# Correlate with average beta values
## Add beta values to data frame
avg_centrality$WB_novelty <- WB_novelty_values
avg_centrality$WB_familiarity <- WB_familiarity_values
## Calculate correlations
r1 <- round(cor(avg_centrality$score1, avg_centrality$WB_novelty), 2)
BF10_1 <- reportBF(correlationBF(avg_centrality$score1, avg_centrality$WB_novelty))
p1 <- pValue(cor.test(avg_centrality$score1, avg_centrality$WB_novelty)$p.value)
r2 <- round(cor(avg_centrality$score1, avg_centrality$WB_familiarity), 2)
BF10_2 <- reportBF(correlationBF(avg_centrality$score1, avg_centrality$WB_familiarity))
p2 <- pValue(cor.test(avg_centrality$score1, avg_centrality$WB_familiarity)$p.value)
scatter1 <- ggplot(avg_centrality, aes(x = score1, y = WB_novelty)) +
geom_point() +
geom_smooth(method = "lm", formula = "y ~ x", colour = base_col[1]) +
#scale_x_continuous(breaks = c(0.50, 0.575, 0.65)) +
scale_y_continuous(breaks = c(-400, -200, 0, 200)) +
coord_cartesian(xlim = c(0.50, 0.65), ylim = c(-400, 200),
expand = FALSE) +
theme_journal() +
labs(x = "Centrality score", y = "Whole-brain novelty\ncontrast (a.u.)")
scatter2 <- ggplot(avg_centrality, aes(x = score1, y = WB_familiarity)) +
geom_point() +
geom_smooth(method = "lm", formula = "y ~ x", colour = base_col[1]) +
scale_y_continuous(breaks = c(-100, 0, 100, 200)) +
coord_cartesian(xlim = c(0.50, 0.65), ylim = c(-100, 210),
expand = FALSE) +
theme_journal() +
labs(x = "Centrality score", y = "Whole-brain familiarity\ncontrast (a.u.)")
##############################################################################_
p <- plot_grid(avg_centrality_plot, scatter1, scatter2,
nrow = 1, labels = "auto", label_size = 10, align = "hv")
# Combine plots
ggsave(p,
filename = paste0(figurePath, "Fig_S05_individual_differences_centrality.png"),
dpi = dpi, width = 180, height = 55, units = "mm")- The correlation between the centrality score and the novelty contrast was r = 0.07, BF10 = 0.35, p = .584.
- The correlation between the centrality score and the familiarity contrast was r = 0.1, BF10 = 0.4, p = .442.
# Load the p-value maps that have been corrected across main effect and correlation
SpaNov1_centrality_p1_file <- "/media/alex/work/Seafile/imaging_results/SpaNov/OLMe_7T_SpaNov_gradient_6lvl_smo4_MSMAll_centrality/cope7.feat/stats/vwc/results_lvl2cope1_dat_ztstat_cfdrp_c3.dscalar.nii"
SpaNov1_centrality_p2_file <- "/media/alex/work/Seafile/imaging_results/SpaNov/OLMe_7T_SpaNov_gradient_6lvl_smo4_MSMAll_centrality/cope7.feat/stats/vwc/results_lvl2cope1_dat_ztstat_cfdrp_c4.dscalar.nii"
# Load the maps
SpaNov1_centrality_p1_xii <- read_cifti(SpaNov1_centrality_p1_file, brainstructures = "all")
SpaNov1_centrality_p2_xii <- read_cifti(SpaNov1_centrality_p2_file, brainstructures = "all")
# Correct again only across the two maps
new_correction <- correct_4_multiple_contrast_xiftis(list(SpaNov1_centrality_p1_xii, SpaNov1_centrality_p2_xii))
# Write to new cifti files
new_name <- str_replace(SpaNov1_centrality_p1_file, pattern = ".dscalar.nii", replacement = "new.dscalar.nii")
write_cifti(new_correction[[1]], cifti_fname = new_name)## Writing left cortex.
## Writing right cortex.
## Writing subcortical data and labels.
## Creating CIFTI file from separated components.new_name <- str_replace(SpaNov1_centrality_p2_file, pattern = ".dscalar.nii", replacement = "new.dscalar.nii")
write_cifti(new_correction[[2]], cifti_fname = new_name)## Writing left cortex.
## Writing right cortex.
## Writing subcortical data and labels.
## Creating CIFTI file from separated components.# Create data frame to compare differences
values <- c(as.matrix(SpaNov1_centrality_p1_xii), as.matrix(SpaNov1_centrality_p2_xii),
as.matrix(new_correction[[1]]), as.matrix(new_correction[[2]]))
compare_p_dist <- data.frame(dist_type = rep(c("Correcting across main effect", "Only correcting across\npositive and negative"), each = length(values)/2),
log_p_values = values)
# Convert from log to actual p-values -log(0.05, 10)
compare_p_dist$p_values <- 10^-(compare_p_dist$log_p_values)
# Get lowest p-value
min_p <- min(compare_p_dist$p_values[compare_p_dist$dist_type == "Only correcting across\npositive and negative"])
# Five lowest p-values
low_5_p <- head(sort(compare_p_dist$p_values[compare_p_dist$dist_type == "Only correcting across\npositive and negative"]))The five lowest p-values are 0.025789, 0.6326812, 0.8176687, 1, 1, 1.
Hippocampal gradient based on maximum statistic
# Create supplemental figure for maximum stat
## Creating a new data frame for combined plotting
HC_data_plotting2 <- HC_data_plotting
HC_data_agg_pos2$Hemisphere2 <- "Average"
HC_supp_fig_data <- data.frame(position = c(HC_data_plotting2$position, HC_data_agg_pos2$position),
max = c(HC_data_plotting2$max, HC_data_agg_pos2$max),
Hemisphere2 = c(HC_data_plotting2$Hemisphere2, HC_data_agg_pos2$Hemisphere2))
## Re-order Hemisphere2 so the order is left, right, then average
HC_supp_fig_data$Hemisphere2 <- factor(HC_supp_fig_data$Hemisphere2,
levels = c("Left hippocampus", "Right hippocampus",
"Average"),
ordered = TRUE)
## Create plot
scatter_max <- ggplot(HC_supp_fig_data, aes(x = -position, y = max)) +
facet_grid(~Hemisphere2, scales = "free_x") +
geom_point(aes(colour = max)) +
geom_smooth(method = "lm", formula = y ~ x, colour = "black") +
labs(x = "Distance to most anterior part (mm)", y = "Novelty Preference") +
scale_x_continuous(breaks = seq(from = -45, to = 0, by = 15), labels = c("45", "30", "15", "0")) +
scale_y_continuous(breaks = 1:6, labels = paste("Lvl", 1:6), limits = c(0.5, 6.5)) +
scale_colour_viridis_c(option = "H", limits = c(1, 6), direction = -1) +
theme_journal() +
theme(strip.background = element_rect(color="white", fill="white")) +
theme(legend.position = "none")
## Save with as 80 mm x 80 mm
ggsave(scatter_max, filename = paste0(figurePath, "Fig_S06_max_scatter.png"),
dpi = dpi, width = 120, height = 40, units = "mm")Hippocampal gradient without averaging
# Create new copy of HC_data and create check sum
HC_data_no_avg <- HC_data
runCodeAgain3 <- check_if_md5_hash_changed(HC_data_no_avg, hash_table_name = "SpaNov_md5_hash_table.csv")
# Run/load the code/data
if(runCodeAgain3){
# Seed
set.seed(21312)
# Other input
lm_formula <- "val ~ position * Hemisphere + tSNR + GS"
nIter <- 100000
colName <- "min"
imageName <- "SpaNov_permut_HC_grad_analysis1_cue-delay_no-avg.RData"
grad_min_permut1_no_avg <- permutation_analysis(HC_data, lm_formula, nIter,
colName, imageName)
# Other input
lm_formula <- "val ~ position + tSNR + GS"
nIter <- 100000
colName <- "min"
imageName <- "SpaNov_permut_HC_grad_analysis1_cue-delay_L_no-avg.RData"
grad_min_permut1_L_no_avg <- permutation_analysis(HC_data[HC_data$Hemisphere == "left", ],
lm_formula, nIter, colName, imageName)
# Other input
lm_formula <- "val ~ position + tSNR + GS"
nIter <- 100000
colName <- "min"
imageName <- "SpaNov_permut_HC_grad_analysis1_cue-delay_R_no-avg.RData"
grad_min_permut1_R_no_avg <- permutation_analysis(HC_data[HC_data$Hemisphere == 'right', ],
lm_formula, nIter, colName, imageName)
} else {
load("intermediate_data/SpaNov_permut_HC_grad_analysis1_cue-delay_no-avg.RData")
grad_min_permut1_no_avg <- results
load("intermediate_data/SpaNov_permut_HC_grad_analysis1_cue-delay_L_no-avg.RData")
grad_min_permut1_L_no_avg <- results
load("intermediate_data/SpaNov_permut_HC_grad_analysis1_cue-delay_R_no-avg.RData")
grad_min_permut1_R_no_avg <- results
}
# Calculate p-value
dist <- grad_min_permut1_no_avg$permuted_values[,5]
critVal <- grad_min_permut1_no_avg$lm$coefficients[6]
grad_min_permut1_no_avg_p5 <- pValue_from_nullDist(critVal, dist, "two.sided")
# Calculate p-value
dist <- grad_min_permut1_L_no_avg$permuted_values[,1]
critVal <- grad_min_permut1_L_no_avg$lm$coefficients[2]
grad_min_permut1_L_no_avg_p1 <- pValue_from_nullDist(critVal, dist, "two.sided")
# Calculate p-value
dist <- grad_min_permut1_R_no_avg$permuted_values[,1]
critVal <- grad_min_permut1_R_no_avg$lm$coefficients[2]
grad_min_permut1_R_no_avg_p1 <- pValue_from_nullDist(critVal, dist, "two.sided")
# Calculate effect sizes
grad_min_permut1_no_avg_eff <- eta_squared(car::Anova(grad_min_permut1_no_avg$lm, type = 2))
grad_min_permut1_L_no_avg_eff <- eta_squared(car::Anova(grad_min_permut1_L_no_avg$lm, type = 2))
grad_min_permut1_R_no_avg_eff <- eta_squared(car::Anova(grad_min_permut1_R_no_avg$lm, type = 2))
# Create data frame
## Create variables
tmp_formula <- c("min ~ position * Hemisphere + tSNR + GS", "min (left) ~ position + tSNR + GS",
"min (right) ~ position + tSNR + GS")
tmp_coef_name <- c("interaction", "position", "position")
tmp_coef_val <- c(grad_min_permut1_no_avg$lm$coefficients[6], grad_min_permut1_L_no_avg$lm$coefficients[2],
grad_min_permut1_R_no_avg$lm$coefficients[2])
tmp_coef_es <- c(grad_min_permut1_no_avg_eff$Eta2_partial[5], grad_min_permut1_L_no_avg_eff$Eta2_partial[1],
grad_min_permut1_R_no_avg_eff$Eta2_partial[1])
tmp_coef_p <- c(grad_min_permut1_no_avg_p5, grad_min_permut1_L_no_avg_p1,
grad_min_permut1_R_no_avg_p1)
## Convert to numeric values and flip
tmp_coef_val <- -as.numeric(tmp_coef_val)
## Add to one data frame
tab_df <- data.frame(Formula = tmp_formula, Coefficient = tmp_coef_name,
beta = tmp_coef_val, eta_squared = tmp_coef_es, p = tmp_coef_p)
## Round columns
tab_df$beta <- round(tab_df$beta, 3)
tab_df$eta_squared <- round(tab_df$eta_squared, 3)
## Create p-values by looping over all values
for(i in seq_along(tab_df$p)){
tab_df$p[i] <-paste("p", pValue(as.numeric(tab_df$p[i])))
}
# Show table
kable(tab_df)| Formula | Coefficient | beta | eta_squared | p |
|---|---|---|---|---|
| min ~ position * Hemisphere + tSNR + GS | interaction | 0.016 | 0.004 | p = .013 |
| min (left) ~ position + tSNR + GS | position | 0.007 | 0.003 | p = .158 |
| min (right) ~ position + tSNR + GS | position | 0.026 | 0.032 | p < .001 |
# Add better names for the hemispheres
HC_data$Hemisphere2 <- ifelse(HC_data$Hemisphere == "right",
"Right hippocampus", "Left hippocampus")
# Create plots
scatter_min_no_avg <- ggplot(HC_data, aes(x = -position, y = min)) +
facet_grid(~Hemisphere2, scales = "free_x") +
geom_jitter(aes(colour = min), height = 0.1, width = 0) +
geom_smooth(method = "lm", formula = y ~ x, colour = "black") +
labs(x = "Distance to most anterior part (mm)", y = "Novelty Preference") +
scale_x_continuous(breaks = seq(from = -45, to = 0, by = 15), labels = c("45", "30", "15", "0")) +
scale_y_continuous(breaks = 1:6, labels = paste("Lvl", 1:6), limits = c(0.5, 6.5)) +
scale_colour_viridis_c(option = "H", limits = c(1, 6), direction = -1) +
theme_journal() +
theme(strip.background = element_rect(color="white", fill="white")) +
theme(legend.position = "none")
## Save with as 80 mm x 80 mm
ggsave(scatter_min_no_avg, filename = paste0(figurePath, "Fig_S07_scatter_no-avg.png"),
dpi = dpi, width = 80, height = 40, units = "mm")Relationship between spatial memory and successful memory encoding
# Calculate the average placement error for each object & subject.
obj_avg <- ddply(OLM_7T_retrieval, c("subject", "objectName"), summarise,
euclideanDistance = mean(euclideanDistance), N = length(subject),
navTime = mean(navTime))
# Rank the objects within each subject
obj_avg <- ddply(obj_avg, c("subject"), mutate, rank = order(euclideanDistance))
obj_avg$f_rank <- factor(obj_avg$rank)
# Visualise results
p1 <- ggplot(obj_avg, aes(x = objectName, y = rank, fill = objectName)) +
geom_jitter(width = 0.25, height = 0, pch = 21) +
geom_boxplot(outlier.shape = NA, alpha = 0.8) +
stat_summary(geom = "point", fun = "mean", col = 'white', size = 0.8, shape = 24, aes(fill = objectName),
position = position_dodge(width = 0.75),
key_glyph = "rect") +
scale_fill_viridis_d(option = "turbo") +
labs(title = "Rank by objects",
x = "Objects", y = "Within-subject rank") +
theme_journal() +
theme(legend.position = 'none',
plot.title = element_text(hjust = 0.5))
p2 <- ggplot(obj_avg, aes(x = f_rank, y = euclideanDistance, fill = f_rank)) +
geom_jitter(width = 0.25, height = 0, pch = 21) +
geom_boxplot(outlier.shape = NA, alpha = 0.8) +
stat_summary(geom = "point", fun = "mean", col = 'white', size = 0.8, shape = 24, aes(fill = f_rank),
position = position_dodge(width = 0.75),
key_glyph = "rect") +
labs(title = "Placement error by rank",
x = "Rank", y = "Avg. placement error (vm)") +
theme_journal() +
scale_fill_viridis_d() +
theme(legend.position = 'none',
plot.title = element_text(hjust = 0.5))
p3 <- ggplot(obj_avg, aes(x = f_rank, y = navTime, fill = f_rank)) +
geom_jitter(width = 0.25, height = 0, pch = 21) +
geom_boxplot(outlier.shape = NA, alpha = 0.8) +
stat_summary(geom = "point", fun = "mean", col = 'white', size = 0.8, shape = 24, aes(fill = f_rank),
position = position_dodge(width = 0.75),
key_glyph = "rect") +
labs(title = "Navigation time by rank",
x = "Rank", y = "Navigation time (seconds)") +
theme_journal() +
scale_fill_viridis_d() +
theme(legend.position = 'none',
plot.title = element_text(hjust = 0.5))
###############################################################################_
# Create base path
base_path <- "/media/alex/work/Seafile/imaging_results/SpaNov/OLMe_7T_encoding1_cue-delay_smo4_MSMAll/cope7.feat/stats/vwc/"
# Load the beta values for encoding contrast
encoding_betaMap_xii <- read_cifti(paste0(base_path, "Y1.dtseries.nii"), brainstructures = "all")
# Get p-value maps
encoding_FDRp_c1_xii <- read_cifti(paste0(base_path, "results_lvl2cope1_dat_ztstat_cfdrp_c1.dscalar.nii"), brainstructures = "all")
encoding_FDRp_c2_xii <- read_cifti(paste0(base_path, "results_lvl2cope1_dat_ztstat_cfdrp_c2.dscalar.nii"), brainstructures = "all")
# Count significant gray-ordinates
# sum(as.matrix(encoding_FDRp_c1_xii) > cutOff) # = 548
# sum(as.matrix(encoding_FDRp_c2_xii) > cutOff) # = 548
# Calculate average successful encoding contrast in novelty masls
encoding_WB_familiarity <- as.matrix(encoding_betaMap_xii)[as.matrix(WB_familiarity_mask) == 1, ]
encoding_WB_familiarity <- colMeans(encoding_WB_familiarity)
encoding_WB_novelty <- as.matrix(encoding_betaMap_xii)[as.matrix(WB_novelty_mask) == 1, ]
encoding_WB_novelty <- colMeans(encoding_WB_novelty)
# Calculate stats
x <- encoding_WB_familiarity
d1 <- signif(mean(x)/sd(x), digits1)
str1 <- mean_SD_str2(x, report_type, digits1, rounding_type)
BF10_1 <- reportBF(ttestBF(x))
x <- encoding_WB_novelty
d2 <- signif(mean(x)/sd(x), digits1)
str2 <- mean_SD_str2(x, report_type, digits1, rounding_type)
BF10_2 <- reportBF(ttestBF(x))
# Create a data frame for plotting
encoding_WB_df <- data.frame(mask = rep(c("Familiarity", "Novelty"), each = length(encoding_WB_familiarity)),
contrast = c(encoding_WB_familiarity, encoding_WB_novelty))
p4 <- ggplot(encoding_WB_df, aes(x = mask, y = contrast, fill = mask)) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_jitter(width = 0.25, height = 0, pch = 21) +
geom_boxplot(outlier.shape = NA, alpha = 0.8) +
stat_summary(geom = "point", fun = "mean", col = 'white', size = 0.8, shape = 24, aes(fill = mask),
position = position_dodge(width = 0.75),
key_glyph = "rect") +
annotate("text", x = 1, y = 200, label = paste("BF10 =", BF10_1), size = 2) +
annotate("text", x = 2, y = 250, label = paste("BF10 =", BF10_2), size = 2) +
labs(title = "",
x = "Mask", y = "Encoding success (a.u.)") +
theme_journal() +
scale_fill_manual(values = cool_warm_colours) +
theme(legend.position = 'none',
plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.margin = unit(c(0, 1, 0, 22.5), "mm"))
# Assemble the plot
p <- plot_grid(p1, p2, p3, NULL, NULL, p4,
align = "h", nrow = 2, #labels = c("a", "b", "c", "d", "", "e"),
label_size = 10, rel_heights = c(1, 1))
# Save
ggsave(p,
filename = paste0(figurePath, "Fig_S08_successful_encoding.png"),
dpi = dpi,
width = 180,
height = 100,
units = "mm")- Encoding effect in familiarity regions: d = -0.068, M = -2.9 (SD = 43), BF10 = 0.16
- Encoding effect in novelty regions: d = 0.43, M = 31 (SD = 73), BF10 = 12.99
Distribution of principal functional connectivity gradient scores across spatial novelty preference levels in the posterior medial cortex
# Compare the PMC gradient with the first connectivity gradient from Margulies et al. (2016)
## Load PMC gradient
PMC_gradient <- read_cifti("cifti_results/SpaNov_gradient_PMC_min_cue-delay.dlabel.nii")
## Add to Margulies_grad data frame
Margulies_grad$PMC_gradient <- c(PMC_gradient$data$cortex_left, PMC_gradient$data$cortex_right)
## Subset to only vertices inside the PMC gradient
Margulies_grad_PMC <- Margulies_grad[Margulies_grad$PMC_gradient != 0, ]
## Make factor
Margulies_grad_PMC$PMC_gradient <- factor(Margulies_grad_PMC$PMC_gradient, levels = 1:6, ordered = TRUE)
## Remove Level5 because there are only 3 vertices
Margulies_grad_PMC <- Margulies_grad_PMC[Margulies_grad_PMC$PMC_gradient != 5, ]
## Create a figure
Margulies_novelty <- ggplot(Margulies_grad_PMC, aes(x = PMC_gradient, y = Grad_1, fill = PMC_gradient)) +
geom_jitter(height = 0, size = 0.5, alpha = 0.1, width = 0.2) +
geom_boxplot(alpha = 0.8, outlier.shape = NA) +
stat_summary(geom = "point", fun = "mean", col = 'white', shape = 24,
position = position_dodge(width = 0.75), size = 0.8) +
scale_fill_manual(values = novFam_gradient[-5]) +
theme_journal() +
theme(legend.position = "none") +
labs(x = "Spatial novelty-familarity gradient (PMC)", y = TeX(r"( $Margulies'\ et\ al^{10}\ 1^{st}\ gradient$ )"))
# Save
ggsave(Margulies_novelty,
filename = paste0(figurePath, "Fig_S09_PMC_gradient_Margulies.png"),
dpi = dpi,
width = 80,
height = 60,
units = "mm")Hippocampal connectivity patterns correlating with Margulies gradient
# Load results
load("data/ignore_fMRI_version1/grad_FRSC/HC_2_cortex_FRSC_Margulies_gradient.RData")
# Apply fisher z-transformation
HC_2_cortex_Margulies_gradient$z_correlation <- z_transform_fisher(HC_2_cortex_Margulies_gradient$correlation)
# Add MNI y coordinate to it
HC_2_cortex_Margulies_gradient$MNI_y <- rep(rep(HC_data$y, 10), 56)
# Divide into anterior and poster
HC_2_cortex_Margulies_gradient$AP <- ifelse(HC_2_cortex_Margulies_gradient$MNI_y >= -21,
"Anterior", "Posterior")
# Calculate average
HC_2_Margulies_agg1 <- ddply(HC_2_cortex_Margulies_gradient, c("Rnum", "Gradient", "AP"),
summarise, z_correlation = mean(z_correlation))
# Make gradient a factor
HC_2_Margulies_agg1$Gradient <- factor(HC_2_Margulies_agg1$Gradient)
# Calculate differences
HC_2_Margulies_agg2 <- ddply(HC_2_Margulies_agg1, c("Rnum", "Gradient"),
summarise, diff = mean(z_correlation[1] - z_correlation[2]))
# Calculate tests for each gradient
BF10 <- rep(NA, 10)
d <- rep(NA, 10)
for(i in 1:10){
# Calculate one sample difference
diff <- HC_2_Margulies_agg2$diff[HC_2_Margulies_agg2$Gradient == i]
d[i] <- mean(diff)/sd(diff)
BF10[i] <- signif(reportBF(ttestBF(diff)), 2)
}
# Make it a data frame for simpler adding
annot_df <- data.frame(label = paste("BF10 =", BF10),
Gradient = factor(1:10),
y = rep(c(0.09, 0.08), 5))
# Visualise
HC_2_margulies <- ggplot(HC_2_Margulies_agg2, aes(x = Gradient, y = diff, fill = Gradient)) +
geom_hline(yintercept = 0, linetype = "dashed", colour = "darkgrey") +
scale_fill_viridis_d() +
geom_boxplot(outlier.shape = NA) +
stat_summary(geom = "point", fun = "mean", col = 'white', shape = 24,
position = position_dodge(width = 0.75), size = 0.8) +
annotate("text", x = annot_df$Gradient, y = annot_df$y, label = annot_df$label, size = 1.5) +
theme_journal() +
theme(legend.position = "none", plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5)) +
labs(x = "Principal gradient", y = "Anterior - posterior difference\nin correlation")
# Save
ggsave(HC_2_margulies,
filename = paste0(figurePath, "Fig_S10_HC_2_Margulies.png"),
dpi = dpi,
width = 120,
height = 60,
units = "mm")Posterior predictive checks for the novelty score model
Find information here: https://mc-stan.org/rstanarm/reference/pp_check.stanreg.html
# Set extra seed because these plot are relatively variable with regard to the tails
set.seed(20240206)
p1 <- brms::pp_check(m_noveltyScore_run1, ndraws = 100) +
scale_colour_manual(values = c("black", base_col[1]), name = "") +
coord_cartesian(xlim = c(-20, 20), expand = TRUE) +
theme_journal() +
theme(axis.line.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
legend.position = "none")
p2 <- brms::pp_check(m_noveltyScore_run2, ndraws = 100) +
scale_colour_manual(values = c("black", base_col[1]), name = "") +
coord_cartesian(xlim = c(-50, 50), expand = TRUE) +
theme_journal() +
theme(axis.line.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
legend.position = "none")
p3 <- brms::pp_check(m_conn2, ndraws = 100) +
scale_colour_manual(values = c("black", base_col[1]), name = "") +
coord_cartesian(xlim = c(-10, 10), expand = TRUE) +
theme_journal() +
theme(axis.line.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
legend.position = "none")
p4 <- brms::pp_check(m_navTime3, ndraws = 100) +
scale_colour_manual(values = c("black", base_col[1]), name = "") +
coord_cartesian(xlim = c(0, 60), expand = TRUE) +
theme_journal() +
theme(axis.line.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
legend.position = "none")
p6 <- brms::pp_check(m_visits_run1, ndraws = 100) +
scale_colour_manual(values = c("black", base_col[1]), name = "") +
coord_cartesian(xlim = c(0, 20), expand = TRUE) +
theme_journal() +
theme(axis.line.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
legend.position = "none")
p7 <- brms::pp_check(m_visits_run2, ndraws = 100) +
scale_colour_manual(values = c("black", base_col[1]), name = "") +
coord_cartesian(xlim = c(0, 20), expand = TRUE) +
theme_journal() +
theme(axis.line.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
legend.position = "none")
p8 <- brms::pp_check(m_timeSince_run1, ndraws = 100) +
scale_colour_manual(values = c("black", base_col[1]), name = "") +
coord_cartesian(xlim = c(0, 500), expand = TRUE) +
theme_journal() +
theme(axis.line.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
legend.position = "none")
p9 <- brms::pp_check(m_timeSince_run2, ndraws = 100) +
scale_colour_manual(values = c("black", base_col[1]), name = "") +
coord_cartesian(xlim = c(0, 2000), expand = TRUE) +
theme_journal() +
theme(axis.line.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
legend.position = "none")
# Solution to change the linewidth found here: https://stackoverflow.com/questions/55450441/how-to-change-size-from-specific-geom-in-ggplot2
p1$layers[[1]]$aes_params$linewidth <- 0.1
p1$layers[[2]]$aes_params$linewidth <- 0.1
p2$layers[[1]]$aes_params$linewidth <- 0.1
p2$layers[[2]]$aes_params$linewidth <- 0.1
p3$layers[[1]]$aes_params$linewidth <- 0.1
p3$layers[[2]]$aes_params$linewidth <- 0.1
p4$layers[[1]]$aes_params$linewidth <- 0.1
p4$layers[[2]]$aes_params$linewidth <- 0.1
p6$layers[[1]]$aes_params$linewidth <- 0.1
p6$layers[[2]]$aes_params$linewidth <- 0.1
p7$layers[[1]]$aes_params$linewidth <- 0.1
p7$layers[[2]]$aes_params$linewidth <- 0.1
p8$layers[[1]]$aes_params$linewidth <- 0.1
p8$layers[[2]]$aes_params$linewidth <- 0.1
p9$layers[[1]]$aes_params$linewidth <- 0.1
p9$layers[[2]]$aes_params$linewidth <- 0.1
# Combine plots
p_comb <- plot_grid(p1, p2, p3, p4, p6, p7, p8, p9,
labels = "auto", label_size = 10, align = "hv")
ggsave(p_comb,
filename = paste0(figurePath, "Fig_S11_ppc.png"),
dpi = dpi,
width = 90,
height = 90,
units = "mm")Distribuion of starting relative to object locations
# Use one participant for visualisation because they are all the same anyway
start_vis_data <- OLM_7T_encoding[OLM_7T_encoding$subject == unique(OLM_7T_encoding$subject)[1], ]
# Remove the gift box that location was variable
no_gift <- obj_locations[obj_locations$objectName != "gift", ]
# Create the plot
start_locations <- ggplot() +
# Circles
geom_polygon(data = circle1, mapping = aes(x = x, y = z), fill = "lightgrey") +
geom_polygon(data = circle2, mapping = aes(x = x, y = z), fill = "darkgrey") +
geom_point(data = start_vis_data, mapping = aes(x = start_x, y = start_z, colour = targetNames)) +
geom_segment(data = start_vis_data, mapping = aes(x = start_x, y = start_z, xend = object_x, yend = object_z,
colour = targetNames)) +
geom_point(data = no_gift, mapping = aes(x = object_x, y = object_z, colour = targetNames), size = 2) +
scale_colour_viridis_d(option = "turbo") +
# Other settings
coord_equal(xlim = c(-90, 90)) +
theme_void() +
theme(plot.background = element_rect(fill = "white", colour = "white")) +
labs(colour = "") +
theme(legend.text = element_text(size=rel(0.5)),
legend.key.height = unit(0.5,"line"),
legend.key.width = unit(0.1,"line"),
legend.title = element_blank(),
legend.box.margin = margin(0, 0, 0, -10))
ggsave(start_locations,
filename = paste0(figurePath, "Fig_S12_start_locations.png"),
dpi = dpi, width = 80, height = 80, units = "mm")Correlation between novelty measures
###############################################################################_
# Number of visits
# Set range
x1_range <- range(m_visits_run1$data$s_onset_rel)
x1_points <- seq_min_2_max(m_visits_run1$data$s_onset_rel)
x2_range <- range(m_visits_run2$data$s_onset_rel)
x2_points <- seq_min_2_max(m_visits_run2$data$s_onset_rel)
# Get subject-level predictions
pred1_df <- predictions(m_visits_run1, newdata = datagrid(s_onset_rel = x1_points, subject = unique),
by = c("s_onset_rel", "subject"))
pred2_df <- predictions(m_visits_run2, newdata = datagrid(s_onset_rel = x2_points, subject = unique),
by = c("s_onset_rel", "subject"))
pred1_df <- as.data.frame(pred1_df)
pred2_df <- as.data.frame(pred2_df)
# Calculate the minimum for each subject, which is used for colouring
pred1_df <- ddply(pred1_df, c("subject"), mutate, min = min(estimate))
pred2_df <- ddply(pred2_df, c("subject"), mutate, min = min(estimate))
# Create plots
## Run 1
p1_1 <- ggplot(data = pred1_df, aes(x = s_onset_rel, y = estimate, group = subject, colour = min)) +
geom_line() +
scale_color_viridis_c() + theme_journal() +
theme(legend.position = "none",
plot.margin = unit(c(1, 2, 1, 2), "mm")) +
labs(title = "Run 1", x = "Time", y = "Number of visits") +
scale_x_continuous(breaks = x1_range, labels = c("Start", "End")) +
scale_y_continuous(breaks = c(2, 6, 10)) +
coord_cartesian(xlim = x1_range, ylim = c(2, 10), expand = FALSE)
## Run 2
p1_2 <- ggplot(data = pred2_df, aes(x = s_onset_rel, y = estimate, group = subject, colour = min)) +
geom_line() +
scale_color_viridis_c() + theme_journal() +
theme(legend.position = "none",
plot.margin = unit(c(1, 2, 1, 2), "mm")) +
labs(title = "Run 2", x = "Time", y = "Number of visits") +
scale_x_continuous(breaks = x2_range, labels = c("Start", "End")) +
scale_y_continuous(breaks = c(6, 10, 14)) +
coord_cartesian(xlim = x2_range, ylim = c(6, 14), expand = FALSE)
###############################################################################_
# Get original scale for time since
# Load
load("event_tables/images/SpaNov_event_file_gradients.RData")
# Convert list to data frame
data <- rbindlist(subj_list, idcol = "subject")
# Add subjects' R number
data$ppid <- subjIDs_R[data$subject]
# Function to match runStartTime
find_runStartTime <- function(ppid, run){
# Get corresponding to find the run in the trial data
anonKey <- lookupTable$anonKey[lookupTable$Rnum == ppid[1]]
# Use the anonKey & run to get runStartTime
runStartTime <- OLM_7T_trial_results$runStartTime[OLM_7T_trial_results$subject == anonKey &
OLM_7T_trial_results$block_num == run[1]]
return(runStartTime[1])
}
data <- ddply(data, c("subject", "ppid", "run"), mutate,
runStartTime = find_runStartTime(ppid, run))
data$subject <- as.character(data$subject)
# Make time relative to the start of the run. The real onset times will be
# slightly different but this will not matter for this.
data$onset_rel <- data$onset - data$runStartTime
# Add run type
data$runType <- "encoding"
data$runType[data$run == 2 | data$run == 4] <- "retrieval"
# Subset to only encoding
data_sub <- data[data$runType == "encoding", ]
# Change run number to match the description in paper. Run 3 is Encoding Run 2
data_sub$run[data_sub$run == 3] <- 2
# Convert run to factor
data_sub$f_run <- as.factor(data_sub$run)
# Run 1
## Remove NA and get correct run
## Use unique name for data frame because of marginaleffects issues
m_timeSince_run1_data <- na.omit(data_sub[data_sub$run == 1, ])
## Scale x and y
m_timeSince_run1_data$s_onset_rel <- scale_m0_sd0.5(m_timeSince_run1_data$onset_rel)
m_timeSince_run1_data$s_timeSince <- drop(scale(m_timeSince_run1_data$lastVisit))
# Use drop() to deal with marginal effects issues
# Run 2
## Remove NA and get correct run
## Use unique name for data frame because of marginaleffects issues
m_timeSince_run2_data <- na.omit(data_sub[data_sub$run == 2, ])
## Scale x and y
m_timeSince_run2_data$s_onset_rel <- scale_m0_sd0.5(m_timeSince_run2_data$onset_rel)
m_timeSince_run2_data$s_timeSince <- drop(scale(m_timeSince_run2_data$lastVisit))
# Use drop() to deal with marginal effects issues
###############################################################################_
# Time since
# Set range
x1_range <- range(m_timeSince_run1$data$s_onset_rel)
x1_points <- seq_min_2_max(m_timeSince_run1$data$s_onset_rel)
x2_range <- range(m_timeSince_run2$data$s_onset_rel)
x2_points <- seq_min_2_max(m_timeSince_run2$data$s_onset_rel)
# Get subject-level predictions
pred1_df <- predictions(m_timeSince_run1, newdata = datagrid(s_onset_rel = x1_points, subject = unique),
by = c("s_onset_rel", "subject"))
pred2_df <- predictions(m_timeSince_run2, newdata = datagrid(s_onset_rel = x2_points, subject = unique),
by = c("s_onset_rel", "subject"))
pred1_df <- as.data.frame(pred1_df)
pred2_df <- as.data.frame(pred2_df)
# Calculate the minimum for each subject, which is used for colouring
pred1_df <- ddply(pred1_df, c("subject"), mutate, min = min(estimate))
pred2_df <- ddply(pred2_df, c("subject"), mutate, min = min(estimate))
# Create plots
## Run 1
p2_1 <- ggplot(data = pred1_df, aes(x = s_onset_rel, y = estimate, group = subject, colour = min)) +
geom_line() +
scale_color_viridis_c() + theme_journal() +
theme(legend.position = "none",
plot.margin = unit(c(1, 2, 1, 2), "mm")) +
labs(title = "Run 1", x = "Time", y = "Time elapsed (seconds)") +
scale_x_continuous(breaks = x1_range, labels = c("Start", "End")) +
scale_y_continuous(breaks = c(50, 125, 200)) +
coord_cartesian(xlim = x1_range, ylim = c(50, 200), expand = FALSE)
## Run 2
p2_2 <- ggplot(data = pred2_df, aes(x = s_onset_rel, y = estimate, group = subject, colour = min)) +
geom_line() +
scale_color_viridis_c() + theme_journal() +
theme(legend.position = "none",
plot.margin = unit(c(1, 2, 1, 2), "mm")) +
labs(title = "Run 2", x = "Time", y = "Time elapsed (seconds)") +
scale_x_continuous(breaks = x2_range, labels = c("Start", "End")) +
scale_y_continuous(breaks = c(50, 300, 550)) +
coord_cartesian(xlim = x2_range, ylim = c(50, 550), expand = FALSE)
###############################################################################_
# Calculate correlation
temp_corr <- cor.test(event_info_df_encoding$visits, event_info_df_encoding$lastVisit)
plot_str <- paste("r", rValue(temp_corr$estimate), "p", pValue(temp_corr$p.value))
# Create scatter plot
corr_novelty_measures <- ggplot(event_info_df_encoding, aes(x = visits, y = lastVisit)) +
geom_point(alpha = 0.1) +
geom_smooth(method = "lm", formula = "y ~ x", colour = base_col[1]) +
annotate("text", x = I(0.8), y = I(0.8), label = plot_str, size = 2) +
theme_journal() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5)) +
labs(x = "Number of visits", y = "Time elapsed (seconds)",
title = "Relationship between novelty measures")
###############################################################################_
# Combine into one figure
combined_plot <- plot_grid(corr_novelty_measures, plot_grid(p1_1, p1_2, p2_1, p2_2,
ncol = 2, labels = c("b", "", "c", ""),
align = "hv", label_size = 10),
ncol = 2, labels = c("a", ""), label_size = 10)
# Save
ggsave(combined_plot,
filename = paste0(figurePath, "Fig_S13_corr_novelty_measures.png"),
dpi = dpi, width = 190, height = 82, units = "mm")